Idealized simulation study of the relationship of disdrometer sampling statistics to the precision of precipitation rate measurement

Due to the discretized nature of rain, the measurement of a continuous precipitation rate by disdrometers is subject to statistical sampling errors. Here, Monte Carlo simulations are employed to obtain the precision of rain detection and rate as a function of disdrometer collection area and compared with World Meteorological Organization guidelines for a one-minute sample interval and 95% probability. To meet these requirements, simulations suggest that measurements of light rain with rain rates R ≤ 0.50 mm h−1 require a collection area of at least 6 cm × 6 cm, and for R = 1 mm h−1, the minimum collection area 5 is 13 cm × 13 cm. For R = 0.01 mm h−1, a collection area of 2 cm × 2 cm is sufficient to detect a single drop. Simulations are compared with field measurements using a new hotplate device, the Differential Emissivity Imaging Disdrometer. The field results suggest an even larger plate may be required to meet the stated accuracy, likely in part due to non-Poissonian hydrometeor clustering.


Introduction
Ground-based precipitation sensors are commonly used to validate remotely sensed precipitation measurement systems, including satellite (TRMM), WSR-88D radar measurements (Kummerow et al., 2000;Fulton et al., 1998), and numerical weather prediction models (Colle et al., 2005) aimed at hydrology, agriculture, transportation, and recreation applications (WMO, 2018;Estévez et al., 2011;Campbell and Langevin, 1995;Brun et al., 1992). The ability of an automated weather station to detect the presence of light pre-cipitation can be crucial for weather forecasting in remote locations where a human observer is not available to verify the presence of rainfall (Horel et al., 2002;Miller and Barth, 2003). Light rain or drizzle can severely impede road safety (Andrey and Yagar, 1993;Andrey and Mills, 2003;Bergel-Hayat et al., 2013;Theofilatos and Yannis, 2014).
Disdrometers measure particle drop size distributions and provide calculated precipitation rate from the integrated mass flux. Among available disdrometers are the mechanical Joss-Waldvogel (JW) disdrometer (Joss and Waldvogel, 1967), laser or optical sensors such as the OTT Parsivel2 (Tokay et al., 2013;Bartholomew, 2014), and video disdrometers such as the 2DVD (Kruger and Krajewski, 2002;Thurai et al., 2011;Brandes et al., 2007). In instrument development, striking a balance between sampling area and accurate fine-scale precipitation detection is a well-known problem. Among other instrument-specific considerations, disdrometer accuracy depends on the sampling area and time interval. Gultepe (2008) highlights a universal difficulty in accurately measuring very light precipitation. In their study, they compared different co-located precipitation sensors and found large absolute and relative errors (32 %-44 %) for all instruments when precipitation rate R < 0.3 mm h −1 , but the relative errors were approximately 10 % when R > 1.5 mm h −1 . Despite these observations, they found that optical and hotplate disdrometers can more accurately detect light precipitation compared to weighing gauges, despite having comparatively small sampling areas, typically 50 cm 2 , where the VRG101 weighing gauge sampling area is 400 cm 2 . Here, we focus solely on the effect of disdrometer sampling area and time interval on the precision of precipitation rate measure-ment distinct from any other uncertainties associated with the instruments themselves.
The World Meteorological Organization (WMO) recommends that a disdrometer measure precipitation with an output averaging time of 1 min (WMO, 2018). Measurement uncertainty is defined as "the uncertainty of the reported value with respect to the true value and indicates the interval in which the true value lies within a stated probability" specified to be 95 %. For liquid precipitation intensity, the uncertainty requirements for rates of 0.2 to 2 mm h −1 are ±0.1 mm h −1 , and for precipitation rates > 2 mm h −1 they are 5 % (see Table 1). For < 0.2 mm h −1 , the requirement is only detection.
Most commercial instruments do not reach such strict standards. Table 1 provides the reported uncertainties for several precipitation instruments. The instrumentation used by Automated Surface Observing Systems (ASOS) weather stations located at major airports utilizes a heated tipping bucket (HTB) for precipitation accumulation and a light-emitting diode weather identifier (LEDWI) for precipitation type and intensity. The LEDWI uses signal power return to determine the drop size distribution of rain or snow and then classifies the precipitation intensity based on the size distribution (Nadolski, 1998). A significant limitation of the ASOS system is that it cannot discriminate drizzle from light rain, and it qualitatively expresses small amounts as a trace (Wade, 2003;Nadolski, 1998).
Vaisala manufactures a range of optical precipitation sensors that detect and categorize precipitation from the forward scattering of a light beam, including the PWD12, PWD22, PWD52, and FD71P. The PWD12 detects rain, snow, unknown precipitation, drizzle, fog, and haze. The PWD22 has the same resolution and accuracy as the PWD12, but it also detects freezing drizzle, freezing rain, and ice pellets (Vaisala, 2019b). The PWD52 has an increased observation range of 50 km compared to 20 km for the PWD22 (Vaisala, 2018a). All three PWD instruments have a precipitation intensity resolution of 0.05 mm h −1 for a 10 min sampling interval at 10 % uncertainty (Vaisala, 2019b). The PWD22 is included with tactical weather instrumentation intended for US military and aviation operations (TACMET) and reports precipitation type in WMO METAR code format (Vaisala, 2018b). The FD71P has a higher stated sampling frequency and resolution than the PWD instruments of 0.01 mm h −1 with 2.2 % uncertainty in a 5 s measurement cycle (Vaisala, 2019a), although there has yet to be independent scientific evaluation of the device.
Hotplate disdrometers offer an alternative with the advantage of requiring fewer assumptions as mass is inferred from the energy required for evaporation. The Yankee Environmental Systems TPS-3100 determines the liquid water precipitation rate of rain or snow by taking the power difference between upward-and downward-facing hotplates as a measure of the latent heat energy required to evaporate precipitation (Rasmussen et al., 2010). The technology is currently marketed as the Pond Engineering Laboratories K63 Hot-plate Total Precipitation Gauge (Pond Engineering, 2020). The hotplate is 5 in. in diameter, or 126.7 cm 2 , which is equivalent to a square with a width of 11.26 cm. It measures precipitation rate with a resolution of 0.10 ± 0.5 mm h −1 and can detect the onset of light snow within 1 min (Yankee Environmental Systems, 2011;Pond Engineering, 2020).
A newer hotplate disdrometer, yet to be commercialized, is the Differential Emissivity Imaging Disdrometer (DEID) developed at the University of Utah. The DEID measures the mass of individual hydrometeors using a hotplate and a thermal camera, which provide accurate, fine-scale measurements. A larger hotplate sampling area increases the operating cost through higher power consumption. The work here was originally motivated by a desire to minimize DEID power and maximize measurement precision, although the calculations are applicable more generally to other disdrometers such as those described. We employ a Monte Carlo approach (Liu et al., 2012(Liu et al., , 2018Jameson and Kostinski, 1999, 2001a to stochastically generate raindrops based on canonical size distributions aimed at determining the minimum required disdrometer collection area and sampling frequency for precise measurement of precipitation rates between 0.01 and 10 mm h −1 . We consider the precipitation rate uncertainty relative to WMO standards. Inherent precipitation measurement uncertainties associated with the instrument mechanism are not addressed here. Where Joss and Waldvogel (1969) approached the problem analytically by assuming the interarrival times of droplets up to 6 mm in diameter are distributed according to a Poisson distribution, here we approach the problem numerically by employing a Monte Carlo approach. In principle, the results should converge, although the Monte Carlo approach also facilitates the calculation here of the time required to measure the "first drop" in a precipitation event. Joss and Waldvogel (1969) defined a sample size as the product of an area and a sample time. Under the assumption that raindrop size follows an exponential distribution, to measure a precipitation rate of R = 1 mm h −1 to a precision of 10 % within 95 % confidence bounds, the required sample size is 1.5 m 2 s, corresponding to a cross-sectional sampling area of A = 250 cm 2 with a square sampling area width W = 15.8 cm for a nominal 60 s collection interval. The required square width W found in this work for the same parameters is 13 cm.

Differential Emissivity Imaging Disdrometer principle
The DEID obtains the mass of individual precipitation particles by assuming conservation of energy during heat transfer from a square plate to a melting hydrometeor. To determine particle size during evaporation, a thermal camera is directed at a heated aluminum sheet. Since aluminum is a thermal reflector (thermal emissivity ≈ 0.03), whereas water is not ( ≈ 0.96), particles have high brightness tempera- ture and appear as white regions on a low brightness temperature, black background. From the measured cross-sectional surface area, temperature, and evaporation time of each hydrometeor, the effective diameter and volume and the mass of each particle can be calculated with high precision . A piece of polyimide tape with ≈ 0.95 is placed on the side of the sampling area as a reference for the differential emissivity calculation and determination of the camera's pixel resolution. For the study described here, the DEID's aluminum plate had an area of 15.24 cm × 15.24 cm, and the camera pixel resolution was 0.2 mm. Polyimide tape applied to the surface restricted the collection area to A ≈ 7 cm × 5 cm, equivalent to a square width of W = 5.8 cm.
Thermal camera imagery with a sampling frequency between 2 and 60 Hz was used to determine the cross-sectional surface area, temperature, and evaporation time of each hydrometeor . From these parameters, individual hydrometeor mass is calculated from conservation of energy, whereby the heat gained by the hydrometeor is equal to the heat lost by the hotplate when evaporating water through where c p is the specific heat capacity of water at constant pressure, T is the difference in temperature between 0 and time t, m is the mass of the hydrometeor, and L eqv is the equivalent latent heat required for the conversion of the hydrometeor to gas. For liquid precipitation L eqv = L v , where L v is the latent heat of vaporization of water. K is the thermal conductivity of the plate, H is hotplate thickness, A(t) is the area of the water droplet at time t, T p is the temperature of the hotplate, and T w (t) is the temperature of the water at time t. When combining the constants into a single value K d , Eq. (1) simplifies to where K d is determined experimentally. The precipitation rate R DEID is calculated from the total mass of hydrometeors evaporated on the hotplate during a given sample time interval. For each frame of the hotplate captured by the thermal camera, where β = 3.6×10 6 mm s m −1 h −1 , f s is the camera resolution (frame s −1 ), A evap is the total area of water on the sampling area (m 2 ), I mean is the pixel intensity related to the temperature difference between the plate and water through T p −T w (t) ≈ (255−I mean )/256×T p , ρ w is the density of water (1000 kg m −3 ), and A hot is the hotplate area (m 2 ).  Hydrometeor catch inefficiency is a large contributor to precipitation rate measurement error, especially in buckettype precipitation gauges (Pollock et al., 2018). Rasmussen et al. (2010) found that the Yankee hotplate had a catch efficiency of only 50 % with a wind speed of 5 m s −1 and 35 % with a wind speed of 8 m s −1 . The DEID underwent a series of calibration wind tunnel tests to determine the effect of wind on mass measurements. During these experiments, mass measurements remained approximately constant, revealing that catch inefficiency is not a contributor to the precipitation rate measurement . The high catch efficiency from the DEID was demonstrated during a storm that took place at Alta, UT, on 16 April 2020 between 00:00 and 16:00 UTC (Fig. 10 of Singh et al., 2021). The colocated weighing gauge was located inside of a wind fence, while the DEID was not. Despite wind speeds during the storm ranging from 4 to 13 m s −1 sustained with 8-19 m s −1 gusts, DEID precipitation measurements were within 6 % of the co-located weighing gauge. For this reason, precipitation rate as a function of catch efficiency is not explored in this work.

Size distribution generation
For a range of collection areas A and time intervals t, a raindrop distribution is generated stochastically and compared with the assumed rate R. Initially, we adopt the Marshall-Palmer (Marshall and Palmer, 1948) drop size distribution where n 0 = 8000 m −3 mm −1 and = 4.1R −0.21 mm −1 . D ranges between 0 and D max in linear bins evenly spaced by D. We establish D max = 6 mm to account for the expected breakup of large raindrops (Villermaux and Bossa, 2009) and D = 0.25 mm for 50 bins. The value of D is arbitrary but was chosen to approximate the spatial measurement resolution of the prototype DEID. For an assumed value of R, the array of drops is stochastically generated according to Eq. (4), where N MP = n(D) D is the total number of drops generated from each bin. Each drop in the bin is assigned a diameter of D mean = D + D/2. For each value of where the coefficient a and the exponent b are determined for the Stokes regime for D mean ≤ 0.08 mm, the intermediate regime for 0.08 mm ≤ D mean < 1.2 mm, and the turbulent regime for D mean ≥ 1.2 mm following Lamb and Verlinde (2011). The maximum value of v is used to determine the sample volume of the generated Marshall-Palmer distribution of drops v max A t m 3 . The calculated size distribution of drops incident on the collection area during sampling time t is then N(D) coll = n(D)Av t, with units of mm −1 . The total calculated number of drops N coll incident on the collection area is a summation of Eq. (6) over the range 0 to D max .

Calculated precipitation rate
N coll drops are randomly sampled assuming the Marshall-Palmer distribution. From the drops that impact the collection area, the calculated precipitation rate R calc is (cf. Lane et al., 2009) where α = 3.6 ×10 −3 m 2 s mm −2 h −1 and N i is the number of drops with diameter D i , with i corresponding to bin number; 100 simulations of R calc are performed for each value of equivalent sampling area width W = √ A, each evenly spaced by 1 cm. The 95th and 5th percentile bounds of R calc are specified as the upper and lower bounds of sampling uncertainty (R calc /R − 1).

First and hundredth drops
To determine the sampling time required to detect the onset of precipitation, 100 simulations were performed for R = 0.01, 0.1, and 1 mm h −1 and for 100 evenly spaced width bins W between 1 and 20 cm. For each width bin, the number size distribution was calculated from n(D) for a 1 m 3 volume directly above the collection area with height h = 1/A. A sample of drops is generated from Eq. (6). To ensure that the drop with the smallest size D mean could feasibly fall from the top to the bottom of the sample volume, t = h/v is maximized using the fall speed of the minimum value of D mean . In general, small drops contribute negligibly to calculations of precipitation rate (Smith et al., 2003), but due to their higher concentrations they may nonetheless be the first detected. Accordingly, the functional form of n(D) is adjusted from the simpler exponential form described by Eq. (4) to a gamma distribution So that small particles with D < 1 mm are not overrepresented (Ulbrich and Atlas, 1984), the drop size distribution is modified by the shape parameter µ and generated according to Eq. (6). Each drop is assigned a random height z above the collection area within a distance h above the plate, and the time elapsed for the plate to detect a drop is t p = z/v. The shortest of these times is the first drop detection time t 1 .
Following the collection approach taken by Marshall et al. (1947), reproduction of the Marshall-Palmer size distribution is assumed to require collection of 100 drops. The time elapsed for the calculated incidence of 100 drops is t 100 . If fewer than 100 drops were obtained in N coll , a new sample of drops is obtained from Eq. (6) with an increased value of t.

Monte Carlo simulations
The sampling uncertainty in the precipitation rate is illustrated in Fig. 2. Precipitation rates of R = 0.02, 0.20, and 2.00 mm h −1 were analyzed for a standard sampling time of 60 s. Collection areas smaller than approximately 6 cm × 6 cm meet WMO standards for R ≤ 0.50 mm h −1 , but a collection area of over 10 cm × 10 cm is required for R > 1 mm h −1 .
To illustrate the relationship between accuracy and sampling time, Fig. 3 compares R = 0.10, 1.00, and 10.00 mm h −1 with sampling times of 10 s, 5 min, and 10 min. For heavy rain rates of 10 mm h −1 , for a time interval of 5 min, a disdrometer sampling area width of 4 cm yields measured rain rates with a precision of ±5 % error Figure 5. Recalculation of R calc (Eq. 7) from 1000 uniformly distributed, randomly sampled iterations of a 1 min DEID dataset from 8 March 2020 at 14:43 MST with rain rate R DEID = 4.2 mm h −1 (a). R calc was recalculated using uniformly distributed random segments of time (c) and hotplate sampling area width (d). As t → 60 and W → W DEID , R calc approaches R DEID .
for 95 % of the measurements. The intersection between 95th and 5th percentile bounds and WMO accuracy criteria occurs in larger collection areas as R increases.
First drop simulation results are shown in Fig. 4. Three simulations were performed using three values of µ, where µ = 0 represents the Marshall-Palmer exponential distribution and closely represents the distributions generated by the uncertainty simulations. Following Ulbrich and Atlas (1984), DEID measurements show values of µ between 1 and 2 best represent the distribution of drops that arrive on the hotplate (Figs. 5b and 6b). Monte Carlo calculations of the size distribution N coll are shown in Fig. 4. Note that the size distribution for N (D) coll does not converge to N (D) MP for long sampling times because smaller drops fall slowly. Rather, the distribution more closely resembles a gamma distribution (Ulbrich and Atlas, 1984). Nevertheless, the contribution of small drops to ground-based measurements of precipitation tends to be small as they are comparatively less massive. Also, precipitation particles form primarily from droplet collisions in the updrafts within clouds and so must attain the size that they fall sufficiently quickly to leave cloud base and fall to the ground where they can be sampled by ground-based instruments (Garrett, 2019). For µ = 2, few of the smallest drops with D < 1 mm are incident on the collection area. A square collection sampling area width of 2 cm × 2 cm is sufficient to detect the onset of light precipitation with a rate of R = 0.01 mm h −1 within 1 min.

Application to DEID measurements
During a field campaign that took place at the University of Utah between April 2019 and March 2020, the DEID recorded the mass and density of individual hydrometeors, along with the 1 min-averaged precipitation rate R DEID for six rain events and five snow events with data spanning 1185 total minutes. DEID particle mass distributions for two contrasting 1 min samples of convective moderate (R DEID = 4.2 mm h −1 ) and light (R DEID = 0.1 mm h −1 ) rain are shown in Figs. 5 and 6. These samples were selected from rain-only events with measured R DEID values closest to the values of R that were analyzed in Sect. 4.1. These disdrometer data are used here in place of an assumed size distribution to calculate R calc (Eq. 7) in 100 iterations, each taken from data randomly sampled over time interval segments of Figure 6. Recalculation of R calc from 1000 iterations of a 1 min DEID dataset from 8 March 2020 at 15:14 MST with rain rate R DEID = 0.1 mm h −1 (a). R calc was calculated using uniformly distributed random segments of time (c) and hotplate sampling area width (d).
As t → 60 and W → W DEID , R calc approaches R DEID . t = 10, 20, 30, 40, and 50 s (Figs. 5c and 6c) and plate area segments of W = 1, 2, 3, 4, and 5 cm (Figs. 5d and 6d). Segments of t include all particles within the DEID collection area, and segments in W encompass the entire 60 s time interval. In a three-dimensional space of plate cross-sectional area and sampling time, the associated number concentration of drops for each case is shown in Figs. 5a-b and 6a-b.
The version of the DEID used in this study had a maximum collection area width of W = 5.8 cm. For moderate rain where R = 1 mm h −1 , the derived minimum required sampling area width to meet WMO requirements is 13 cm. That is, the size of plate used was insufficient for the measurement of rain this intense. For light rain, the statistical uncertainty bounds at 95 % confidence converge to within ±0.1 mm h −1 of the DEID 1 min measured rate for a plate collection area width of W ∼ 3 cm. This value is larger than the 2 cm suggested by the Monte Carlo calculation shown in Table 2. One possibility is that the raindrop interarrival time and spacing were not in fact Poissonian (Jameson and Kostinski, 2001b). In the event of clustering, a larger plate would be required to provide an accurate assessment of the average rain rate during a given time interval.
To assess whether the presence of non-Poissonian clustering is the case, two-point correlation functions η were calculated following Shaw et al. (2002). A value of unity indicates interarrival times that are Poissonian and values greater than unity the presence of non-random clustering. Based on the location of hydrometeor centroids as they arrive on the DEID plate, storm-averaged values of η were found to be equal to 1.01 for rain falling under light winds on 8 March 2020, between 13:38 and 15:49 MST, and equal to 1.10 for the period between 05:00 and 07:34 MST earlier that day when high winds were present. For contrast, the value of η was 1.55 for a snow event with large aggregate snowflakes that took place on 14 January 2020 at 12:43-14:06 MST. While more extensive analysis is required, the implication is that non-Poissonian clustering can occur.

Conclusions
A Monte Carlo approach was used to determine the minimum required cross-sectional collection area for a disdrometer to measure a given precipitation rate with a WMO target precision at 95 % probability for a 1 min collection period. Intrinsic instrument uncertainties were not considered, only those associated with statistical sampling errors associated with the raindrop size distribution. Following these criteria, a square collection area of 6 cm × 6 cm is sufficient to detect the onset of light rain with R ≤ 0.50 mm h −1 .
For R > 1 mm h −1 , a sample area of over 10 cm × 10 cm is required, although a smaller collection area may achieve the required accuracy by increasing the sampling time. For example, in 10 min, a 4 cm × 4 cm collection area can measure 10 mm h −1 precipitation rates to within the WMO required precision 95 % of the time. A collection area as small as 2 cm × 2 cm may detect the onset of light drizzle with R = 0.01 mm h −1 within 1 min, even in instances where small particles in the drop size distribution fall too slowly to intercept the collection area.
Theoretical results obtained from Monte Carlo simulations were compared with observed field measurements from a new precipitation sensor, the Differential Emissivity Imaging Disdrometer, for both light and moderate rain. Randomly selected segments of decreasing sampling time and area from the DEID were used to recalculate the precipitation rate. The results suggest a larger plate may be required to meet a specified precision than those indicated by the Monte Carlo simulations that were performed. A possible explanation is the presence of non-Poissonian clustering that was revealed by two-point correlation function calculations, particularly during high wind and snow events. The results presented here have general implications for the sampling limitations of other widely used particle-by-particle disdrometers such as the PARSIVEL with a sample area of 48.6 cm 2 or effective W of 7 cm (Battaglia et al., 2010) and the 2DVD (Kruger and Krajewski, 2002) with a W of 10 cm. Despite their sizable collection areas, like the Differential Emissivity Imaging Disdrometer, they may nonetheless fail to meet WMO standards if operated at a nominal 1 min sampling interval. Specific heat capacity of water at constant pressure = 4.28 × 10 3 J K −1 kg −1 D max Maximum generated drop diameter = 6 mm D mean Mean diameter value in each bin (mm) f s Camera resolution (frame s −1 ) H Hotplate thickness = 0.0508 m K Thermal conductivity of Al = 205 W m −1 K −1 K d Experimentally derived constant = 1.54 ×10 −3 kg s −1 K −1 m −2 L eqv Equivalent latent heat required for the conversion of the hydrometeor to gas L v Latent heat of vaporization of water = 2.26 ×10 6 J kg −1 W DEID Square sampling area width of the DEID (cm) α Conversion factor = 3.6 × 10 −3 m 2 s mm −2 h −1 β Conversion factor from m s −1 to mm h −1 = 3.6 × 10 6 mm s m −1 h −1 Thermal emissivity Slope parameter = 4.1R −0.21 for the Marshall-Palmer distribution ρ w Density of liquid water = 1000 kg m −3