The objective of this study is to characterize the signature of dynamical microphysical processes on reflectivity–rainfall (Z–R) relationships used for radar rainfall estimation. For this purpose, a bin model with explicit microphysics was used to perform a sensitivity analysis of the shape parameters of the drop size distribution (DSD) as a function of time and rainfall regime. Simulations show that coalescence is the dominant microphysical process for low to moderate rain intensity regimes (R < 20 mm h−1) and that the rain rate in this regime is strongly dependent on the spectral properties of the DSD (i.e., the shape). The time to equilibrium for light rainfall is at least twice as long as in the case of heavy rainfall (1 h for stratiform vis-à-vis 30 min for thunderstorms). For high-intensity rainfall (R > 20 mm h−1), collision–breakup dynamics dominate the evolution of the raindrop spectra. The time-dependent Z–R relationships produced by the model converge to a universal Z–R relationship for heavy intensity rainfall (A = 1257; b ∼ 1) centered on the region of Z–R space defined by the ensemble of over 100 empirical Z–R relationships. Given the intrinsically transient nature of the DSD for light rainfall, it is proposed that the vertical raindrop spectra and corresponding rain rates should be modeled explicitly by a microphysical model. A demonstration using a multicolumn simulation of a Tropical Rainfall Measuring Mission (TRMM) precipitation radar (PR) overpass over Darwin for a stratiform event during the Tropical Warm Pool–International Cloud Experiment (TWP-ICE) is presented.
In radar meteorology, Z–R relationships are used to estimate rainfall intensity from reflectivity measurements. The Z–R relationships relate the value of the measured reflectivity (Z) to the value of the rain rate (R) according to the general formula (Marshall et al. 1947)
where the reflectivity (Z) is expressed in millimeters to the sixth power per meter cubed and the rain rate (R) is in millimeters per hour. Empirical estimates of the two power law coefficients (the prefactor A and the exponent b) vary according to meteorological, climatic, and physiographic conditions. The most common Z–R relationships worldwide were compiled by Battan (1973) for the period 1947–70 and subsequently by Raghavan (2003) for Z–R relationships published after 1970. Overall, about 100 Z–R relationships are available in the literature for different types of rain events (stratiform rainfall, convective rainfall, showers, thunderstorms, monsoon, tropical storms, hurricanes, etc), climatic regions (tropics, higher latitudes, specific location, etc), and season of the year. Attempts to organize Z–R relationships based on rainfall characteristics (Fujiwara 1965; Ulbrich 1983) have met only partial success by identifying an apparent coherence among orographic and thunderstorm precipitation with regard to the location of the power law coefficients in the (A, b) parameter space (Uijlenhoet and Pomeroy 2001). The large variability of the empirically determined power-law coefficients highlights the intrinsically local (transient) nature of the relationships, both in space and time, as storm systems evolve and propagate. An additional source of variability (and uncertainty) stems from the use of different sensors (rain gauges, radar, disdrometers, etc.) and rain estimation techniques; different Z–R relationships can be obtained for the same event using different types of measurements and sampling procedures (Stout and Mueller 1968; Smith et al. 1993; Ciach and Krajewski 1999; Amitai 2000; Campos and Zawadzki 2000; Steiner and Smith 2000; Tokay et al. 2001). The effect of scale difference between radars and rain gauges is another relevant factor of uncertainty, as pointed out by Steiner and Smith (2004). This raises the specter of calibration to specific observing systems, which further undermines the generality of empirical Z–R relationships. Nevertheless, the simplicity of Eq. (1) is appealing and it continues to motivate research toward deriving general Z–R relationships.
Steiner et al. (2004) examined the dependency of Z–R relationships on the shape of time-invariant analytical expressions of the drop size distribution (DSD) (e.g., exponential, gamma, monodisperse). In addition to weather conditions, geographical location, and measurement processes, one factor that fundamentally impacts the estimation of Z–R relationships is the transient evolution of rainfall spectral properties due to dynamical microphysical processes (coalescence, breakup, accretion, evaporation, and condensation). In this work, we focus on transients and address microphysical processes explicitly. Specifically, the research goal is to investigate the link between small-scale physical processes such as coalescence and breakup that occur at hydrometeor scale and Z–R relationships used for rainfall intensity estimation over the characteristic measurement scale of radars. For this purpose, a microphysical model for the evolution of the drop spectra in the presence of coalescence and breakup (Prat and Barros 2007a,b, hereinafter PB07a and PB07b) was used to generate thousands of raindrop spectra for a wide range of rainfall regimes using different DSDs as boundary and initial conditions. A parametric study was conducted to access the influence of microphysical processes on time-varying DSD properties, including the time necessary to achieve equilibrium under steady boundary conditions and the statistics of drop–drop interactions.
The paper is organized as follows: First, a review of the microphysical model, as well as basic formulas and interpretations of DSD and empirical Z–R relationships, is provided in section 2. Second, the analysis of model simulations is presented in section 3. Section 4 examines the relationship between transient DSD properties and microphysical processes in the light of drop–drop interactions statistics conditional on the rain rate regime. An explicit microphysical simulation of vertical profiles of reflectivity and rain rate for a Tropical Rainfall Measuring Mission (TRMM) precipitation radar (PR) (http://trmm.gsfc.nasa.gov) overpass is presented in section 5. The manuscript wraps up with a summary and conclusions in section 6.
2. Microphysical and statistical models
a. Microphysical model and integral properties
In this work, we use a dynamical microphysics bin model in two configurations: 1) a zero-dimensional (0D) single box (PB07a) and 2) a 1D column (rain shaft) (PB07b). Briefly, the bin model uses a conservative scheme (number and mass) for the discretization of the general continuous stochastic collection–breakup equation with state-of-the-art aggregation and breakup kernels (PB07a; PB07b; Prat et al. 2008, hereinafter PBW08). The discretization of the raindrop spectra was implemented on an optimized irregular grid (PB07b) with 40 bin categories (drop size intervals) with diameters ranging from 0.01 to 0.7 cm. Initial and boundary conditions were specified by discretizing a specified DSD into the 40 bin categories above. In the case of the 1D implementation of the model, boundary conditions were applied at the top of the rain shaft. For the applications described in this manuscript, the boundary conditions and time step (Δt = 1 s) are maintained constant during each simulation.
The integral properties of the simulated raindrop spectra [i.e., the rain rate (R) and reflectivity (Z)] are computed as follows:
where Ni(z, t) is the drop number density in each bin category expressed by unit of diameter (cm−4). The fall velocity Vi (cm s−1) for a drop of diameter di is calculated according to Best (1950). Finally, the reflectivity (dBZ) is obtained by
b. Statistical models—Drop size distributions
Three types of formulations are used typically to express the raindrop size spectra: exponential (Marshall and Palmer 1948, hereinafter MP48), gamma (Ulbrich 1983), and normalized gamma (NG) distributions (Willis 1984; Testud et al. 2001; Bringi et al. 2004). The expression for the exponential raindrop size distribution N(D) (cm−4) is given by (MP48)
where N0 (cm−4) and Λ (cm−1) are the parameters of the DSD. Based on their observations, Marshall and Palmer further proposed a constant value N0 = 0.08 cm−4 and the following relationship between the slope Λ and the nominal rain rate RMP48 (mm h−1):
Ulbrich (1983) proposed the gamma distribution:
where N0 is expressed in cm−4−μ.
Finally, the use of a normalized gamma distribution to describe the shape of the DSD has gained increased popularity (Willis 1984; Testud et al. 2001; Illingworth and Blackman 2002; Bringi et al. 2004):
where NW is the normalized intercept parameter (cm−4), D0 is the median volume diameter (cm), f (μ) is the shape factor function, and μ is the shape factor:
Hereinafter, we focus on the more general form of the normalized gamma distribution for analysis. The ranges selected for the gamma distribution parameters (μ, NW, D0) are similar to those found elsewhere (Gorgucci et al. 2006):
In addition, the rain rate is constrained to
c. Empirical Z–R relationships and the DSD
For the empirical Z–R relationships reported in the literature (Battan 1973; Raghavan 2003), parameters A and b vary respectively between 16.6 < A < 865 and between 1.08 < b < 2.87, with most of the expressions (75%) falling in the range 150 < A < 550 and 1.19 < b < 1.7. Ulbrich (1983) proposed a classification of the power law coefficients in the (A, b) space according to rainfall type (stratiform, orographic, showers, thunderstorms). Uijlenhoet and Pomeroy (2001) noted that there is large variability among empirical Z–R relationships for isolated showers and orographic and thunderstorm rainfall, whereas Z–R relationships for stratiform events display the highest degree of correlation. We address this point in section 3 by mapping the regions of DSD parameter space generated by the microphysical model that are physically consistent with a number of selected Z–R relationships.
When looking at the rain rate (R) and radar reflectivity factor (Z) in the space (NW, D0) for two extreme values of the shape factor (μ = −1, μ = 5), the shape factor (μ) appears to have little to no influence on the calculation of the rain rate [Eq. (2)] for low to moderate values (R < 50 mm h−1). However, there is a significant difference in the calculation of the reflectivity [Eq. (3)] corresponding to different values of the shape parameters. They are clearly distinct from each other for reflectivities of 20 dBZ and higher, which is close to reflectivity thresholds of 18–20 dBZ commonly used in operational radar algorithms as a rain/no rain flag. This asymmetry in the Z–R relationship for light rainfall may be one source of uncertainty (ambiguity) in radar retrieval algorithms. For the remainder of this manuscript, the parameter space (μ, NW, D0) is limited to rainfall regimes in the 0.1 mm h−1 < R (mm h−1) < 500 mm h−1 range.
3. Bulk rainfall properties
The modeling results presented in sections 3a and 3b were obtained with the single box model (PB07a). The goal was to investigate the link between microphysical processes in warm rain and empirical Z–R relationships for a wide range of initial conditions and rainfall regimes. Thousands of initial raindrop spectra [NG or Marshall–Palmer (MP) distributions] were generated, and the corresponding bulk rain rate and reflectivity were calculated at each time step. The modeling results presented in section 3c were obtained with the 1D column (rain shaft) model (PB07b), and the analysis is conducted for raindrop spectra and Z–R relationships at the bottom of the rain shaft, corresponding to surface rainfall conditions.
a. Matching transient DSDs and empirical Z–R relationships
Figure 1 displays the normalized rain rate (Rr = Rt=equ./Rt=0; left) and the normalized reflectivity (Zr = Zt=equ./Zt=0; right), where the normalization is with respect to the initial conditions for two values of the shape factor (μ = −1 and μ = 5) at equilibrium. Steady state is achieved when relative variations of the rain rate and for two moments of the DSD [the drop number concentration (zeroth moment) and the reflectivity (sixth moment with respect to the diameter)] are less than 0.01% (0.0001) between two consecutive time steps (Δt = 1 s). The liquid water content, which is proportional to the third moment with respect to the diameter, remains constant throughout the simulation because of global mass conservation constraints in the model and the fact that evaporation was switched off. Visual inspection of the graphs indicates that the median volume diameter D0 is the DSD parameter that represents coalescence efficiency (Rr > 1; Rr < 1 for breakup), whereas there is little to no sensitivity to NW in the transient evolution of rain rate (horizontal lines) except for light rainfall (D0 < 0.1 cm). Similar conclusions can be drawn for reflectivity. For shape factor μ ≤ 2 (not shown), Zr increases with μ for the same D0, which indicates that the number concentration of small drops is directly associated with the shape factor. Concurrent increases in N0 and μ suggest there is predominance of breakup in the low μ range.
Figure 2 summarizes the time necessary to achieve steady state using the box model. The results are presented by using an NG DSD with μ = 0 as the initial condition that is discretized into the 40 class categories of the bin model. Note that the isolines for the time to equilibrium are closely aligned with the rain rate isolines and that the initial shape of the DSD (NW, D0, μ) has only a limited effect on the time to equilibrium (for equilibrium time less than 1 h). The shape factor plays little role in the time necessary to achieve equilibrium and the results obtained for the two extreme values of the shape factor (μ = −1 versus μ = 5) are of the same order of magnitude (not shown), with a shorter time to reach equilibrium in the case of an initial DSD with a high number of drops (μ = −1 versus μ = 5). Note that for very light (stratiform) rainfall (R < 10 mm h−1), steady state is reached only after 1 h of simulation, which highlights the transient nature of short-lived rain systems as compared to the characteristic time scales of measurement.
To explore the transient nature of the DSD, the time-varying DSDs simulated by the model as well as related integral properties are contrasted with commonly used empirical relationships. The basic idea is to identify the parameters (NW, D0, μ) of the initial NG DSD [Eq. (11)] that evolve under the influence of coalescence and breakup in the model and at some point in time a specific Z–R relationship [Eq. (1)]. In other words, the objective is to map the regions in parameter space that correspond to the ensemble of initial NG DSDs that lead to DSDs that match known (desired) empirical Z–R relationships after some simulation time. This exercise required over 7000 model simulations. The regions of parameter space (NW, D0, −1 < μ < 5) were mapped for four different types of rainfall as in Ulbrich (1983): stratiform (Fig. 3a), thunderstorm (Fig. 3b), orographic (Fig. 3c), and showers (Fig. 3d). Regardless of μ, and for all rain types with exception of just one case for orographic rainfall, the region of parameter space when R > 50 mm h−1 is only sparsely populated and is never reached when NW and D0 are both high. Indeed, higher rainfall rates (R > 50 mm h−1) are only achieved for NW < 2 × 104 (mm m3)−1 and D0 > 0.15 cm, indicating a relatively sparse DSD population of very large drops typical of tropical rain. For most of the cases presented here, the intersection of regions of parameter space (NW, D0, −1 < μ < 5) occupied by different Z–R relationships (i.e., the regions where any Z–R relationship will do), and thus where determination of DSD parameters is ambiguous, corresponds to rain rate intensity of R < 20 mm h−1 for stratiform precipitation and thunderstorms (Figs. 3a,b), R < 50 mm h−1 for orographic precipitation (Fig. 3c), and R < 5 mm h−1 for showers (Fig. 3d). These results can provide guidance with regard to the need to specify physically meaningful bounds to DSD parameters for different Z–R relationships. The empirical Z–R relationship for the stratiform case (Z = 313R1.25) spans almost all the parameter space (NW, D0, −1 < μ < 5) with the exception of a small domain located around NW < 2.5 × 104 (mm m3)−1 and D0 ≈ 0.2 cm. This behavior is observed for values of the exponent b ≈ 1.25 regardless of the value of the prefactor A. The region occupied by the empirical Z–R relationship (Z = 31R1.71) proposed by Blanchard (1953) corresponds to higher values of the rain rate R > 50 mm h−1 (and up to 500 mm h−1), which are not typical of orographic rainfall unless there is embedded convection, perhaps localized enhancement of rainfall from a band associated with a mesoscale vortex (e.g., Lang and Barros 2002) or a tropical cyclone, for example. Another explanation is feeder–seeder interactions between stacked orographic clouds leading to the enhancement of low-level rainfall (Barros and Lettenmaier 1993). The feeder–seeder mechanism is similar from a microphysics perspective to extraordinarily effective coalescence in a quickly evolving storm environment.
Interestingly, the time necessary to verify a specific Z–R relationship with the bin model is quasi-instantaneous (inferior to 2 min) when the initial condition provided is the DSD corresponding to Blanchard’s Z–R relationship (A = 31), increases to slightly less than 30 min for thunderstorms, and varies between 30 min to 1 h for stratiform rain, which is consistent with the early results provided in Fig. 2 using rain rate as an indicator of rainfall regime. The time necessary to reach any specific Z–R relationship decreases with decreasing values of μ because a larger population of small droplets tends to increase the efficiency of collision–breakup processes and thus lengthen the duration of the transient period for increasing μ. These differences in the characteristic time scale of the dynamics of microphysical processes for each rain event suggest that to derive useful empirical Z–R relationships from observations, the latter should be integrated over time scales corresponding to the rain rate regime at the time of observation. Nevertheless, doing so will not remove all uncertainties because of the heterogeneity of rainfall in space and time and the challenge in transferring information between the operational scales of radars and rain gauges.
b. Transient Z–R relationships
A sense of the transient behavior of Z–R relationships can be acquired from model-simulated DSDs obtained by imposing as initial condition a Marshall–Palmer DSD for different values of the nominal rain rate (RMP48 = 1, 2, 5, 10, 20, 50, 100, 200, and 300 mm h−1; Fig. 4a). Note that after specifying the initial conditions using the MP DSD, points in the Z–R diagram are plotted using the actual value of the rain rate computed over all the bins according to Eq. (2), and not the nominal rain rate RMP48 as in Eq. (3b). Otherwise, the initial parameters for A and b at t = 0 would follow the MP48 Z–R relationship for stratiform rain (Z = 200R1.6). For reference, all Z–R relationships in the published literature are also displayed in Fig. 4a as background (gray lines). The time-evolving Z–R relationship for the MP DSD sweeps different regions and in different directions the common Z–R space for rain rates above and below ∼20–30 mm h−1.
Specifically, the evolution of the DSD toward equilibrium depends on the value of the initial nominal rain rate of the MP DSD: for rain rates R < 30 mm h−1, the evolution is toward increasing values of rain rate and reflectivity as the DSD evolves under the combined influence of drop coalescence and breakup mechanisms, and the opposite is true for higher values. In particular, the larger drops in the DSD spectra vanish by breakup and a higher number of small drops is created; conversely, a higher number of drops is created by coalescence in the midrange drop diameter (0.02–0.03 cm) domain, which plays an important role in the value of the DSD integral properties such as the rain rate (the velocity-weighted third-order moment with respect to the drop diameter) and reflectivity (proportional to the sixth-order moment with respect to the drop diameter; PB07a). This behavior changes for rain rates above 30 mm h−1, veering in the opposite direction toward decreasing values of rain rate and reflectivity.
The Z–R parameters corresponding to different stages of this transient behavior of the MP DSD were obtained via a simple linear fit in log–log space and are presented in Table 1. Under the combined effects of coalescence and breakup, and for increasing simulation time toward steady state, the parameter A increases while b decreases. Previously, Wilson and Brandes (1979) reported increasing A and decreasing b for a case of coalescence only, and decreasing A and b for cases of aerodynamic breakup only.
When an ensemble of NG DSDs is used to specify initial conditions (Fig. 4b), we find that the Z–R relationships simulated by the bin model collapses to a single line for rain rates above 30 mm h−1, exhibiting the same behavior as discussed for the MP DSD in the context of Fig. 4a. Finally, the time to equilibrium decreases with increasing rain rate as expected (PB07a). Thus, for heavier rainfall (R > 20 mm h−1), the numerical experiments suggest that a regime characterized by a unique Z–R relationship is attained quickly, independent of initial and boundary conditions, whereas at lower rates ambiguity remains for longer periods of time and thus so does the corresponding uncertainty in radar rainfall retrieval algorithms. The model-based Z–R relationship for high-intensity rainfall (A = 1257; b ∼ 1) is centered on the region of Z–R space defined by the ensemble of over 100 empirical Z–R relationships for the same rainfall intensity regime (R > 30 mm h−1). Note that the parameters are close to the values obtained at steady state independent of the rain rate regime (Table 1) as expected: heavier rainfall reaches steady state faster if all else remains equal. This result confirms List (1988), who found a linear Z–R relationship (b = 1) in the case of steady tropical rain. It also confirms that once the effects of coalescence and breakup are balanced, rain rate and reflectivity are directly proportional, regardless of the coalescence and breakup kernels used (see PB07a). Uijlenhoet et al. (2003) also reported experimental evidence of a linear Z–R relationship (b = 1) in the case of rainfall events with very high rain rates (R > 100 mm h−1) and comparable prefactor values (A = 1750; A = 1540; A = 853). Jameson and Kostinski (2001a,b) associated linear Z–R relationships (b = 1) with statistically homogeneous rain and argued that power law Z–R relationships (b ≠ 1) were not justified physically and could be attributed to errors in the observations and statistically inhomogeneous data. In the very light rainfall regime (R < 1 mm h−1), model simulations suggest much higher actual reflectivity than that implied by the empirical Z–R relationships for the same rain rate. This is likely to be a measurement bias at low rain rates.
This study highlights the need to introduce rainfall intensity thresholds to cluster and interpret observations from field experiments in order to capture the relevant microphysics. By simply fitting all observations with one Z–R relationship independently of rain rate regime, slow and fast microphysics are mixed, which would explain the very different model-based (A, b) parameters as compared to the literature.
In operational rainfall retrieval algorithms, optimum parameters are typically determined by minimizing a cost function that relates reflectivity and rain rate. However, it is critical that physically meaningful constraints reflecting microphysical behavior be considered to obtain realistic parameters (μ, NW, D0). Results obtained with the single box model provide valuable input concerning the limits for the domain of variability of parameters that are estimated to retrieve the shape of the DSD in radar rainfall estimation algorithms. Based on the characteristics of the rain event observed (stratiform, orographic, thunderstorm), comprehensive boundaries can be set on the DSD parameters in order to retrieve the actual DSD.
c. Transient evolution in the Z–R space for a one-dimensional shaft model
Using static boundary conditions imposed at an altitude of 4000 m in a vertical one-dimensional rain shaft, DSD evolution throughout the shaft was simulated as the drops are free to settle by gravity alone (no vertical wind speed is specified); that is, the rain shaft is isolated from its environment. At each time step, and at each height, the DSD is modified by coalescence and breakup mechanisms. The evolution of the DSD and related integral properties (R, Z) at the bottom of the shaft after the profiles reach steady state are examined next. The underlying assumption is that we are dealing with DSDs corresponding to stable raindrops that can only change by interaction with other equally stable raindrops (coalescence and breakup). As showed by Carbone and Nelson (1978), the DSD spectra can experience significant changes during the lifetime of a storm system depending on moisture convergence and storm dynamics. Recently, PBW08 showed that by using transient boundary conditions (DSDs retrieved from vertically pointing radar), the model was able to simulate with high fidelity the temporal evolution of a stratiform rainfall event as compared to radar and ground-based disdrometer measurements. Here, the objective is to characterize the influence of the shape of initial DSD on DSD evolution when undergoing changes due to coalescence and breakup alone. That is, we focus the evolution of the DSD within a rain shaft that is isolated from the storm environment.
Simulations were performed (not shown) to account for the influence of μ (=1, 0, 2, 5) for rain rates ranging for 0.7 to 300 mm h−1. As pointed out previously (Fig. 4), the value of the shape factor μ is more important at low to moderate rain intensity (R ≤ 10 mm h−1) and negligible at higher rain rates (R > 30 mm h−1). Relative to previous results obtained for different parameter combinations (μ = 0, NW, D0), the shape factor μ is less sensitive to the dynamics of the microphysical processes than the other DSD parameters (NW, D0). Smaller standard deviations are observed for different values of μ than for different values of NW and D0 regardless of the value of the rain rate. Table 2 summarizes the value of the coefficients (A, b) of the Z = ARb relationship at the bottom of the rain shaft at the end of 1-h simulation as a function of the exponent (μ) (thus, not at equilibrium) and independently of rain rate. The parameter A decreases with increasing value of μ, whereas the exponent b increases with increasing value of μ as in stratiform rain conditions (Ulbrich 1983). This is consistent with the fact that the homogenous 1D model as used here is comparable to a realistic situation with no, or limited, updrafts/downdrafts, as is the case in stratiform conditions. For simulation of realistic thunderstorm-like conditions to investigate the universal Z–R relationship, a different model setup including vertical winds would be necessary.
Figure 5 displays the results for the 1D model in Z–R space for different values of NW and D0 and for μ = 0. Simulations were performed for 45 cases representing nine different rain rate intensities (0.5, 1, 2, 5, 10, 20, 50, 100, 200 mm h−1). For each one of the nine selected nominal rain rates, the value of NW is fixed (NW = 4000, 10 000, 20 000, 40 000, 60 000 mm−1 m−3) and D0 is adjusted accordingly (ΔD0 = 0.005 cm), as reported in Table 3. All simulations are run for a total time of 3600 s. Starting from uniform initial conditions at each level throughout the column (same DSD at all z levels), this duration is long enough to ensure that equilibrium profiles of the integral properties (drop number concentration, liquid water content, rain rate, and reflectivity) are achieved. It is important to make a distinction between the time necessary to reach equilibrium in the case of a 0D model (for which the DSD at the previous time step is injected at the next time step and so on) and the case of a 1D model in which the same DSD is imposed at the top of the column at each time step. In the first case, we obtain an equilibrium DSD; in the second, we obtain equilibrium profiles.
Three empirical Z–R relationships [stratiform (MP48), orographic (Blanchard 1953), and thunderstorm (Jones 1956)] are plotted for perspective. The standard deviation computed at ground level indicates that for similar rain rates, the shape of the initial DSD (μ = 0, NW, D0), plays an important role at low to moderate rain rates as indicated by the dispersion of results in the Z–R diagram for R < 20 mm hr−1. For higher rain rates (R ≈ 30 mm h−1) the simulation results collapse to one line in the Z–R diagram. Figure 5 therefore confirms the split in behavior noted in section 3 with regard to the temporal evolution of the Z–R relationship as a function of rainfall intensity, although the rain rate threshold for the regime with universal Z–R relationship appears to have shifted toward higher rainfall rates.
Overall, simulation results obtained with the single box and shaft models are consistent with observations reported elsewhere (Atlas and Chmela 1957; Atlas 1964; Tokay and Short 1996; Atlas et al. 1999) that showed a reduction in the dispersion of measured Z–R for rain rates between 20 and 50 mm h−1. In addition, for R > 30 mm h−1, we obtain a comparable linear Z–R relationship (b = 1.06) than with the single box model but with a lower prefactor value (A = 844) closer to the value (A = 742) reported by List (1988). As mentioned earlier, those differences can be explained by the different collisional breakup kernel formulations used in both studies as well as by the different methods used for the resolution of the general stochastic collection–breakup equation (PB07a).
Although the Z–R relationship for stratiform rainfall agrees well with the mean model Z–R relationship for R < 10 mm h−1, it fails to capture the simulated model behavior for higher rainfall rates. Again, these results suggest that instead of using storm type as a means to classify Z–R relationships, the microphysics warrant the use of rainfall intensity thresholds to cluster and analyze observations.
To account for the dynamics of microphysical processes throughout the shaft, simulation results (not shown) indicate that the standard deviation with respect to the vertical dimension of the shaft [i.e., computed at both ends of the shaft; between ground level and the top (h = 4000 m) where the boundary conditions (DSDs) are imposed] is almost constant regardless of the value of the rain rate. The value of the standard deviation with respect to the vertical dimension is comparable with the standard deviation obtained with respect to the shape of the DSD (NW, D0) for a rain rate intensity approximately R ≤ 10 mm h−1. This further confirms the predominance of the influence of the spectral properties of the DSD in explaining the transient evolution of the Z–R relationship for low to moderate rain rates (R ≤ 20 mm h−1). On the contrary, the evolution of the DSD is closely tied to rain rate for higher values of the rain rate (R ≥ 20 mm h−1). At higher rain rates, an increased number of drop–drop interactions leading to coalescence and breakup should lead to a faster evolution of the drop spectra. This is examined next.
4. Transient DSD characteristics and drop–drop collisions statistics
To understand the role of drop interactions and microphysical processes in determining whether a specific Z–R relationship is controlled by the initial shape of the DSD (light rainfall) or the rain rate (heavy rainfall), we now focus on the statistics of drop–drop collisions. One key question concerns the threshold frequency of coalescence–breakup interactions necessary to be in one or another regime. The raindrop collision rate (CTΔt) between two drops of size DL and DS with clustering effects is given by Eq. (15) in McFarquhar (2004, hereafter MF04):
where VL and VS are the fall velocities (Atlas et al. 1973) of drops DL and DS, respectively, and Ω(t) is the drop pair cross-correlation function given by MF04. In the first part of the paper, the expression from Best (1950) was used for the drop fall velocity. Although there is no effective difference between model simulations obtained with fall velocities calculated using either the Best (1950) or Atlas et al. (1973) equations (see, e.g., Guzel and Barros 2001), we use Atlas et al. in this section because it can be used to derive directly analytical expressions for the computation of the raindrop collision rate (CTΔt):
with Ω0 = 1.3 and m = 0.055 s−1.
The raindrop collision rate expresses the total number of raindrop collisions occurring within a parcel of the atmospheric domain by unit of time. Here the raindrop collision rate is Δt = 1 s, computed for the generalized gamma DSD as a function of (μ, NW, D0) along with results for the particular case of a MP DSD (μ = 0, NW = 8000 mm−1 m−3, D0 = 3.67/Λ) is shown in Fig. 6. The total number of drop–drop collisions (CT) increases with increasing values of NW and D0 and decreasing values of the shape factor μ. Moreover, regardless of the rain rate intensity, CT increases with the increasing size of the small drop population (i.e., decreasing μ and D0), and regardless of NW, differences in CT observed for different values of the shape factor (μ = −1, 0, 2, 5) increase with increasing D0 and rain rate intensity. Interestingly, a MP DSD with varying values of the nominal rain rate captures the mean behavior in terms of total drop collision CT when compared to the entire set of normalized gamma DSDs that were considered in this study with parameters defined by Eqs. (11) and (12).
For a MP DSD, the total number of collisions CT is expected to increase with increasing rain rate intensity R (McFarquhar and List 1991; MF04), whereas in the case of a more general normalized gamma DSD, our results indicate that the shape factor μ that controls the shape of the DSD in the small drops domain (but D0 and NW influence it as well because those parameters are not totally independent) has a major influence on the value of CT. That is, the shape (μ, D0, NW) of the DSD is a critical factor determining the nature of drop–drop interactions and collision–breakup dynamics.
Finally, simulation results obtained with both the box model and the shaft model have shown that a rain rate value of about 20 mm h−1 separates two microphysical regimes. Above this value, steady state can be achieved within a reasonable time, and the initial DSD has no influence on the equilibrium situation. For lower rain rates (R < 2 mm h−1), the equilibrium DSD is only achieved in the asymptotic sense for very long durations and microphysical processes are controlled by the parameters of the DSD (NW, D0, μ). Despite the regime transition between DSD-controlled (for lower rain rate) and rain rate–controlled (for higher rain rate) regimes in Z–R space, this analysis confirms that rain intensity is a surrogate index of an increase in the drop–drop interaction probability. Moreover, we note that the limit R ≥ 20 mm h−1 determined previously, and which marks the boundary between DSD shape and rain rate controls, corresponds to a total number of raindrop collisions (CT) in the range 1 < CT (m−3 s−1) < 100.
The transition between DSD-controlled and rain rate–controlled regimes observed here by using a microphysical model with detailed representation of coalescence and breakup mechanisms is consistent with results reported elsewhere (Uijlenhoet et al. 2003; Steiner et al. 2004). Steiner et al. (2004) provided analytical Z–R relationships as a function of the form of the DSD (exponential, gamma, monodisperse) and distinguished three microphysical conditions referred to as number controlled, size controlled, and a mix of both conditions. The results obtained here provide quantitative information regarding the boundary of those regimes: the rain-rate (or number)-controlled and shape-controlled regimes are located respectively above and below the value of R = 20 mm h−1.
5. Application of the microphysical model to radar rainfall retrieval
The 1D model (PB07b; PBW08) was extended to multicolumn geometry using space–time-varying TRMM PR reflectivity products (1C21 and 2A25) as a top boundary condition to estimate the surface rain rate at the time of the satellite overpass. In particular, we present here results for the same stratiform event (23 January 2006) simulated by PBW08 using vertically pointing radar (VPR) reflectivity at 4-km altitude from the Tropical Warm Pool–International Cloud Experiment (TWP-ICE) as the top boundary condition. Reflectivity from 1C21 (Zm) corresponds to the effective radar reflectivity factor at 13.8 GHz without any propagation loss (correction due to rain or any other atmospheric gas). The TRMM PR algorithm 1C21 performs a conversion of the power and noise value to reflectivity (Zm) in case of rain without rain attenuation correction. Reflectivity from 2A25 (Zcorr) includes a correction for attenuation loss. The TRMM PR algorithm 2A25 employs a hybrid of the surface reference method and the Hitschfeld–Bordan method (Iguchi et al. 2000) and provides a rain rate estimate for each beam. The differences between Zm and Zcorr correspond therefore to the attenuation correction and increase with reflectivity level (Fig. 7). The measured reflectivity field (Zm) used here is for the same day (23 January 2006) and for the time of the TRMM satellite overpass (around 2310 UTC) at the VPR site (12.44°S, 130.96°E).
An area of 100-km diameter around the VPR site was selected for analysis, which corresponds to 407 TRMM vertical profiles and thus 407 model columns (rain shafts). For each column, the TRMM time-varying reflectivity at 4 km (Fig. 8a) was converted into a Marshall–Palmer DSD corresponding to stratiform conditions as indicated in the TRMM PR product 2A23 that provides information regarding rain type (classified as stratiform, convective, or other) for this particular event (not shown). Figure 8a shows the reflectivity field used to specify initial boundary conditions at the top of the 4-km shafts for this exercise. The surface rain rate field (h = 0 km) produced by the model for a 15-min simulation is presented in Fig. 8b. The selection of a 15-min simulation represents a compromise between the difficulty of emulating an instantaneous situation with a snapshot of the reflectivity field and longer simulation times (1 h) that did not present significant differences. For intermediate to light rainfall rates (R < 10 mm h−1), the reflectivity ratio (Rh=4km/Rh=0km) between the top level (h = 4 km) and the ground (h = 0 km) (not shown) is below one (0.92 < Rh=4km/Rh=0km < 0.99), consistent with the predominance of coalescence processes in DSD evolution dynamics as noted earlier in the context of Fig. 5.
Using the measured reflectivity (Zm) from 1C21 as a boundary condition, the average residuals computed for the 407 TRMM vertical profiles (included in an area of D = 100 km around the VPR site) between the model and TRMM PR 2A25 (not shown) are 3%–6% for the reflectivity expressed in reflectivity decibel scale at different heights in the rain shaft and 11% for the rain rate (maximum of 150%) at ground level (Fig. 9a). The rain rate is calculated in the TRMM 2A25 algorithm by integrating an NG DSD that reproduces the reflectivity (Zcorr). The spread in the relative error for R < 5 mm h−1 is expected based on our discussion in section 3 regarding the ambiguity of Z–R relationships for light rainfall. Overall, the residual average values are of the same order of magnitude as those reported by PBW08, who presented an intercomparison of numerical results obtained with the 1D model, VPR measurements at several heights inside the 1D shaft, and surface measurements performed with a Joss–Waldvogel disdrometer (JWD). However, a detailed analysis of the spatial patterns of the TRMM PR rain rate (2A25) versus model rain rate at ground level shows a different behavior depending on the profile location, with the largest differences observed for profiles located above the ocean as compared with profiles located over land and/or coastline (Fig. 9b). Specifically, for profiles located over land or coastline, the modeled rain rate at ground level is about 20% higher than the TRMM 2A25 rain rate; by contrast, the modeled rain rate is about 40% lower than the TRMM 2A25 rain rate for profiles over ocean and increases for rain rate R > 7–10 mm h−1. These discrepancies may be explained, at least in part, by the radar operating at 13.8 GHz or nonuniform beam-filling corrections in the TRMM 2A25 algorithm (Iguchi et al. 2000). Briefly, in the 2A25 algorithm the path-integrated attenuation (PIA) is performed in two ways using the Hitschfeld–Bordan (HB) technique and the surface reference technique (SRT) (Meneghini et al. 2000) and uses a maximum likelihood estimate to determine the optimum PIA (Iguchi et al. 2000). For example, there appears to be a systematic difference in the effective attenuation correction between land and ocean (Fig. 9c), which might be explained by the approach used to calculate PIA. The selection of a priori parameters to compute the attenuation correction factor or the differences observed between the two techniques (HB and SRT) at low to moderate rain rates as mentioned elsewhere (Robertson et al. 2003) may introduce discrepancies. This point, however, is beyond the scope of the present study.
Nevertheless, a simple intercomparison of the surface rain rate estimated by the model using Zm (1C21) and Zcorr (2A25) as top boundary conditions suggests that there may indeed be significant differences in the reflectivity corrections applied over the ocean and over land between 1C21 and 2A25 products, which are reflected in the significant differences in the regression lines in Figs. 10a and 10b. In particular, note the threshold rainfall implicit from the intercept in Fig. 10b for the regression analysis over the ocean. The Z–R relationships derived from 2A25 TRMM reflectivity and rain rate are comparable with the MP Z–R relationship (Z = 200R1.6) for profiles over land/coastline (Z = 216R1.44). For the profiles over the ocean, the Z–R relationship obtained by regression between 2A25 TRMM reflectivity and rain rate is Z = 396.44R0.82, with a value for the exponent (b = 0.82) smaller than any other values reported in the literature (Battan 1973; Raghavan 2003). Analysis of the region of parameter space (NW, D0, −1 < μ < 5) that verifies these Z–R relationships over land/coast (Z = 216R1.44) and over ocean (Z = 396R0.82) shows that the latter corresponds to the domain where R < 5–10 mm h−1 (Fig. 11). Furthermore, for rain rates R > 10 mm h−1, no DSD, regardless of the formulation used (NG or MP), fits the ocean Z–R relationship according to the microphysical simulations. Of course, these Z–R relationships were derived in this case study based on the data available and do not necessarily fall in line with the hierarchy of transformations in the operational 2A25 algorithm. Nevertheless, this example points out one potential use of the model to support operational retrieval as a means to establish physically meaningful constraints and eliminate ambiguity in the retrieval of light to moderate rainfall rates.
6. Summary and conclusions
The study consisted of characterizing the signature of dynamical microphysical processes on Z–R (reflectivity–rainfall) relationships used for radar rainfall estimation. A parametric study performed with a microphysical model (PB07a) allowed us to reverse-engineer the DSD characteristics and contributions of different microphysical processes (coalescence and breakup) that verify well-established Z–R relationships for particular configurations (thunderstorms and orographic and stratiform precipitation). Numerical results indicate that only DSDs with rain rates above 10 mm h−1 achieve steady state within 1 h or less. In addition, nonlinear effects were observed for the transient evolution of modeled Z–R relationships for an MP DSD with a change of dynamic for nominal rain rates around 20–30 mm h−1.
Simulations performed with a 1D rain shaft model (PB07b) showed how the DSD evolution throughout the shaft depends on the shape of the initial DSD for low to moderate rain intensity (<10 mm h−1), whereas the shape of the DSD plays little role for higher values of rain rate (>20–30 mm h−1). A series of numerical experiments performed on statistics of drop–drop interactions confirmed the fact that the shape of the DSD was key in the microphysical process dynamics rather than the rain rate intensity for light and moderate rainfall. However, a significant number of drop–drop interactions [total number of raindrop collisions: 1 < CT (m−3 s−1) < 100] is necessary to modify the DSD under the combined influence of coalescence and breakup within an acceptable time duration (<30 min) and moderate to high fall distance (2–4 km). This result indicates that most of the rainfall regimes encountered in nature correspond to transient conditions.
For high-intensity rainfall regimes (R > 20 mm h−1), collision–breakup dynamics dominate the evolution of the raindrop spectra, and the sensitivity of rain rate to the shape of the DSD is not important. Indeed, analysis of the time-dependent Z–R relationships produced by the model suggests convergence to a universal Z–R relationship for heavy-intensity rainfall. The model-based Z–R relationship (A = 1257; b ∼ 1) is centered on the region of Z–R space defined by the ensemble of over 100 empirical Z–R relationships for the same rainfall intensity regime. One implication of these findings is that to capture the relevant microphysics in interpretive studies of field experiments and thus reduce ambiguity for light rainfall, observations should be analyzed according to rainfall intensity rather than rainfall class type.
Consistent with the ambiguity brought to bear by the use of specific Z–R relationships to estimate rainfall in light to moderate rainfall regimes, our analysis of simulations with column model simulation using TRMM reflectivity to specific boundary conditions at the top of the rain shaft shows there is greater variability in model versus algorithm estimates for light rainfall (R < 7 mm h−1). We also found significant differences in rain estimates over land and over the ocean where the TRMM PR rain rate can be up to 80% higher than the model in the domain R > 7–10 mm h−1. This case study, which compares TRMM PR products with results obtained with a multicolumn configuration of the bin microphysical model, illustrates the potential use of the microphysical model to estimate surface rain rate over large areas. Under the condition of availability of reliable reflectivity (or DSD) measurements at a given level of the atmospheric column, the bin microphysical model could be used to obtain the value of the surface rain rate without the use of empirical Z–R relationships or assumptions concerning the shape of the DSD, or, for example, to estimate integral properties where measurement confidence is lower (e.g., below 2 km), where PBW08 have shown the model can do very well as compared to disdrometer data. Nevertheless, there is still a crucial need for fundamental experimental work concerning the microphysical aspect of drop–drop interactions (Barros et al. 2008), especially in the small diameter range of the DSD.
This research was supported in part by NASA Grant NNX07AK40G (for the second author). We thank Dr. Eyal Amitai and Dr. Ali Tokay for useful discussions, and the reviewers for valuable comments and suggestions.
Corresponding author address: Dr. Ana P. Barros, Duke University, Box 90287, 2457 CIEMAS Fitzpatrick Bldg., Durham, NC 27708. Email: email@example.com