Statistical characteristics of raindrop size distribution during rainy seasons in the Beijing urban area and implications for radar rainfall estimation

Raindrop size distribution (DSD) information is fundamental in understanding the precipitation microphysics and quantitative precipitation estimation, especially in complex terrain or urban environments which are known for complicated rainfall mechanism and high spatial and temporal variability. In this study, the DSD characteristics of rainy seasons in the Beijing urban area are extensively investigated using 5-year DSD observations from a Parsivel² disdrometer located at Tsinghua University. The results show that the DSD samples with rain rate < 1 mm h⁻¹ account for more than half of total observations. The mean values of the normalized intercept parameter (log ₁₀N_w) and the mass-weighted mean diameter (D_m) of convective rain are higher than that of stratiform rain, and there is a clear boundary between the two types of rain in terms of the scattergram of log ₁₀N_w versus D_m. The convective rain in Beijing is neither continental nor maritime, owing to the particular location and local topography. As the rainfall intensity increases, the DSD spectra become higher and wider, but they still have peaks around diameter D∼0.5 mm. The midsize drops contribute most towards accumulated rainwater. The D_m and log ₁₀N_w values exhibit a diurnal cycle and an annual cycle. In addition, at the stage characterized by an abrupt rise of urban heat island (UHI) intensity as well as the stage of strong UHI intensity during the day, DSD shows higher D_m values and lower log ₁₀N_w values. The localized radar reflectivity (Z) and rain rate (R) relations (Z=aR^b) show substantial differences compared to the commonly used NEXRAD relationships, and the polarimetric radar algorithms R(K_dp), R(K_dp, Z_DR), and R(Z_H, Z_DR) show greater potential for rainfall estimation.

1 Introduction

Raindrop size distribution (DSD) provides fundamental information on precipitation microphysics. Understanding the DSD variability is of great importance in remote sensing observations of precipitation and microphysical parameterizations in numerical weather prediction (NWP) models. For example, the DSD serves as a fundamental bridge in deriving the Z–R relationships used by ground-based weather radar (Battan, 1973; Uijlenhoet and Stricker, 1999) and spaceborne radar (i.e., TRMM PR: Iguchi et al., 2000; and GPM DPR: Hou et al., 2014) for quantitative precipitation estimation (QPE). The NWP systems coupled with various DSD models can capture more detailed horizontal and/or vertical rainfall information so as to improve the accuracy of precipitation predictions (Abel and Boutle, 2012; Fadnavis et al., 2014; McFarquhar et al., 2015; Saleeby and Cotton, 2004). In addition, the DSD is highly related to the kinetic energy of rainfall that has substantial impact on the soil erosions (Angulo-Martinez and Barros, 2015; Caracciolo et al., 2011; Ellison, 1945; Kinnell, 2005; Lim et al., 2015), which is critical to further understanding of the runoff processes and mitigation of subsequent flood hazards (Angulo-Martinez and Barros, 2015; Smith et al., 2009).

Numerous studies have been devoted to the statistical characteristics of DSD worldwide. It is found that the DSD characteristics vary with geographical locations, climate regimes, seasons, rain types, and even diurnal cycles (Dolan et al., 2018; Ji et al., 2019; Seela et al., 2018; Tokay and Short, 1996; L. Wen et al., 2017). Dolan et al. (2018) classified the global DSD characteristics into six groups by analyzing 12 global disdrometer datasets across three latitudes using principal component analysis. They found that the physical processes shaping the DSD characteristics were likely to vary as a function of location. The comparison of DSD in northern and southern China in Tang et al. (2014) showed that there was a clear difference in precipitation microphysical parameters between different regimes during convective rain, while the difference was less notable for stratiform events. The DSD analysis in Beijing (G. Wen et al., 2017) and Taiwan (Seela et al., 2018) also indicated that there were significant differences in DSD between summer and winter rainfall, and both showed the diurnal variation. In addition, the DSD may exhibit high variability in special weather systems. For example, DSD of the tropical cyclones has a higher concentration of small and middle size drops as well as a lower mass-weighted mean diameter (i.e., D_m) in all types of rain compared with the non-tropical cyclone in Darwin (Deo and Walsh, 2016).

Beijing, the capital of China, is a very densely populated metroplex with a population higher than 21 million. It is more vulnerable to extreme weather events such as torrential rainfall and floods (Zhang et al., 2013). Since the hydrology response in urban areas is sensitive to the spatial and temporal variability of rainfall (Cristiano et al., 2017), rainfall monitoring networks with high temporal and spatial resolution (e.g., dense network of automatic weather stations by de Vos et al., 2017; remote sensing network described by Chen and Chandrasekar, 2015, and Cifelli et al., 2018) have been established in several metropolitan areas. The rapid urbanization and complex topography have further exacerbated the high variability of precipitation in the Beijing urban area, posing challenges to precipitation observations and forecast (Song et al., 2014; Yang et al., 2013a, 2016). This also highlights the importance of understanding local DSD characteristics to better quantify the urban precipitation.

Several studies on DSD characteristics in the Beijing area have been conducted. Tang et al. (2014) studied the DSD characteristics and the polarimetric radar parameters for convective and stratiform rain from July to October 2008 and compared with other regions using a first-generation laser-based optical Particle Size and Velocity (Parsivel¹) disdrometer manufactured by OTT Messtechnik, Germany. G. Wen et al. (2017) investigated the statistical properties of summer and winter precipitation in Beijing, including the bulk properties, raindrop fall velocity, axis ratio, and DSD, using a two-dimensional video disdrometer (2DVD) and a micro-rain radar (MRR). Ji et al. (2019) analyzed the microphysical structure of DSD using 14-month DSD measurements from a second-generation Particle Size and Velocity (Parsivel²) disdrometer in Beijing.

However, these studies are mainly focused on summer time (June–September or July–October) or with very limited measurements from one season or two, which are not sufficient to represent local DSD characteristics, especially the monthly variability, during the rainy seasons ranging from May to October. In addition, the impacts of the urban heat island (UHI) effect on rainfall microphysical properties have never been studied in the literature, as the DSD measurements used in previous studies are more likely collected in the suburban area.

This paper presents a comprehensive study of DSD properties using 5-year (2014–2018) continuous observations in the Beijing urban area, aiming to advance our understanding and characterizations of DSD in urban regions as well as parameterization in remote sensing retrievals and NWP models. The DSD properties, their variabilities, as well as the potential applications in radar QPE are detailed.

This paper is organized as follows. Section 2 describes the dataset and methods for data quality control and analysis. The characteristics of DSD parameters for all rainfall events combined, different rainfall types, different rain rate classes, different periods of a day, and different months are detailed in Sect. 3. Section 4 presents the implications for radar QPE and the parameterization errors of different DSD-based radar rainfall algorithms. Summary and conclusion are given in Sect. 5.

2 Data and method

2.1 Dataset

In this study, a Parsivel² disdrometer is used, which is deployed at Tsinghua University campus, Beijing, China (hereafter referred to as THUD). Figure 1 illustrates the specific location of THUD (40.002^∘ N, 116.324^∘ E; 91 m a.s.l. – above sea level) relative to the Beijing metroplex. It is an optical disdrometer with a 54 cm² horizontal sample area, and it is configured with 1 min sampling resolution to measure the DSD and fall velocity of raindrops (Löffler-Mang and Joss, 2000). The velocity and particle sizes are divided into 32 non-uniform bins, varying from 0.05 to 20.8 m s⁻¹ for velocity and 0.062 to 24.5 mm for particle diameter.

Figure 1 (a) The topography of Beijing and (b) the locations of DSD studies in the Beijing area; the red mark represents the location of the Parsivel² disdrometer deployed at Tsinghua University campus in this study and the green and purple makers represent locations in the studies by G. Wen et al. (2017) and Ji et al. (2019), respectively. The map data are available under the Open Database License. © OpenStreetMap contributors 2019. Distributed under a Creative Commons BY-SA License.

The DSD measurements were collected from June 2014 to December 2018. Lyu et al. (2018) compared the accumulated rainfall computed from the DSD data with the rainfall measurements from an automatic weather station 350 m away from THUD to cross-check the reliability of both instruments. Since most rainfall in the Beijing area occurs during the warm season, which usually lasts from May to the end of October (Song et al., 2014), this study uses the data collected from May to October to analyze the DSD characteristics.

2.2 Method

The direct measurements from the disdrometer are the number of raindrops at each velocity (i) and diameter (j) bin. Here, we take the mid value of each bin as the corresponding value. Then the maximum diameter D_max (mm) of raindrops can be obtained directly from the data and the total number of raindrops T_d can be calculated:

$$\begin{array}{}\text{(1)}& T_{\mathrm{d}}=\sum _{i=\mathrm{1}}^{\mathrm{32}}\sum _{j=\mathrm{1}}^{\mathrm{32}}{n}_{i,j},\end{array}$$

where n_i,j stands for the drop number at each bin.

The number concentration of raindrops per unit volume for the jth diameter bin can be calculated as follows:

$$\begin{array}{}\text{(2)}& N\left(D_j\right)=\sum _{i=\mathrm{1}}^{\mathrm{32}}{\displaystyle \frac{n_{i,j}}{A\cdot \mathrm{\Delta }t\cdot V_i\cdot \mathrm{\Delta }{D}_j}},\end{array}$$

where D_j (mm) is the mid value of the jth diameter bin; N(D_j) is in m⁻³ mm⁻¹; A is the sampling area in m²; Δt is the sampling time interval in s; A and Δt are, respectively, 0.0054 m² and 60 s in this study; ΔD_j (mm) is the diameter spread for the jth diameter bin; V_i (m s⁻¹) is the mid-value fall speed for the ith velocity class.

Because of the measurement error, especially for larger size drops (Tokay et al., 2014), the empirical terminal velocity–diameter (V–D) relationship in Atlas et al. (1973) is adopted in this study:

$$\begin{array}{}\text{(3)}& V\left(D_j\right)=\mathrm{9.65}-\mathrm{10.3}\exp \left(-\mathrm{0.6}{D}_j\right).\end{array}$$

The Gamma model (Ulbrich, 1983) in the following form has been proven to be suitable for describing the raindrop spectra.

$$\begin{array}{}\text{(4)}& N\left(D\right)=N_{\mathrm{0}}{D}^{\mathit{\mu }}\exp (-\mathrm{\Lambda }D),\end{array}$$

where D (mm) is the raindrop diameter; N(D) (mm⁻¹ m⁻³) is the number concentration of raindrops per unit volume per diameter interval; N₀  (mm${}^{-\mathrm{1}-\mathrm{\mu }}$ m⁻³), μ and Λ are the scale, shape and slope parameters of the Gamma distribution, and these three parameters can be derived using gamma moments (GMs) (Kozu and Nakamura, 1991; Tokay and Short, 1996) or maximum likelihood methods (Montopoli et al., 2008). When μ=0, the Gamma form DSD degenerates into an exponential DSD model.

In this study, we use the normalized gamma DSD described by Testud et al. (2000) to describe the natural variations of DSD (Bringi and Chandrasekar 2001; Dolan et al., 2018).

$$\begin{array}{}\text{(5)}& N\left(D\right)=N_{\mathrm{w}}f\left(\mathit{\mu }\right){\left({\displaystyle \frac{D}{D_{\mathrm{m}}}}\right)}^{\mathit{\mu }}\exp \left[-(\mathrm{4}+\mathit{\mu }){\displaystyle \frac{D}{D_{\mathrm{m}}}}\right],\end{array}$$

where N_w (m⁻³ mm⁻¹) is the normalized intercept parameter; D_m (mm) is the mass-weighted mean diameter. N_w, D_m, and f(μ) are calculated as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{m}}={\displaystyle \frac{\sum _{j=\mathrm{1}}^{\mathrm{32}}N\left(D_j\right)\cdot D_j^{\mathrm{4}}\cdot \mathrm{\Delta }{D}_j}{\sum _{j=\mathrm{1}}^{\mathrm{32}}N\left(D_j\right)\cdot D_j^{\mathrm{3}}\cdot \mathrm{\Delta }{D}_j}},\text{(7)}& {\displaystyle }& {\displaystyle }{N}_{\mathrm{w}}={\displaystyle \frac{{\mathrm{4}}^{\mathrm{4}}}{\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}}}\left({\displaystyle \frac{{\mathrm{10}}^{\mathrm{3}}W}{D_{\mathrm{m}}^{\mathrm{4}}}}\right),\text{(8)}& {\displaystyle }& {\displaystyle }f\left(\mathit{\mu }\right)={\displaystyle \frac{\mathrm{6}(\mathrm{4}+\mathit{\mu }{)}^{\mathit{\mu }+\mathrm{4}}}{{\mathrm{4}}^{\mathrm{4}}\mathrm{\Gamma }(\mathit{\mu }+\mathrm{4})}}.\end{array}$$

The integral parameters of total number concentration N_t (m⁻³), rain rate R (mm h⁻¹), liquid water content W (g m⁻³) and the mass spectrum standard deviation σ_m (mm) are also calculated in this study based on the following equations.

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{N}_{\mathrm{t}}=& \int N\left(D\right)dD=\int N_{\mathrm{w}}f\left(\mathit{\mu }\right){\left({\displaystyle \frac{D}{D_{\mathrm{m}}}}\right)}^{\mathit{\mu }}\\ & \exp \left[-(\mathrm{4}+\mathit{\mu }){\displaystyle \frac{D}{D_{\mathrm{m}}}}\right]dD=\sum _{i=i}^{\mathrm{32}}\sum _{j=\mathrm{1}}^{\mathrm{32}}{\displaystyle \frac{n_{i,j}}{A\cdot \mathrm{\Delta }t\cdot V_i}},\end{array}\text{(10)}& {\displaystyle }\begin{array}{rl}R& ={\displaystyle \frac{\mathrm{6}\mathit{\pi }}{{\mathrm{10}}^{\mathrm{4}}{\mathit{\rho }}_{\mathrm{w}}}}\sum _{j=\mathrm{1}}^{\mathrm{32}}V\left(D_j\right)D_j^{\mathrm{3}}N\left(D_j\right)\mathrm{\Delta }{D}_j\\ & =\sum _{j=\mathrm{1}}^{\mathrm{32}}R\left(D_j\right)\mathrm{\Delta }{D}_j,\end{array}\text{(11)}& {\displaystyle }W={\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{6}\times {\mathrm{10}}^{\mathrm{3}}}}\sum _{j=\mathrm{1}}^{\mathrm{32}}{D}_j^{\mathrm{3}}N\left(D_j\right)\mathrm{\Delta }{D}_j,\text{(12)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{m}}=\sqrt{{\displaystyle \frac{\sum _{j=\mathrm{1}}^{\mathrm{32}}{\left(D_j-D_{\mathrm{m}}\right)}^{\mathrm{2}}N\left(D_j\right)\cdot D_j^{\mathrm{3}}\cdot \mathrm{\Delta }{D}_j}{\sum _{j=\mathrm{1}}^{\mathrm{32}}N\left(D_j\right)\cdot D_j^{\mathrm{3}}\cdot \mathrm{\Delta }{D}_j}}},\end{array}$$

where ρ_w is the water density (1.0 g cm⁻³); R(D_j) (mm h⁻¹ mm⁻¹) is the rain rate at the jth diameter class, and it is normalized by the total rain rate R as $R(D_j{)}^{\mathrm{norm}}=\frac{R\left(D_j\right)}{R}$ in the analysis to resolve the contribution of different raindrop sizes to the rainfall intensities. The median volume diameter D₀ (mm) is defined such that drops smaller than D₀ contribute to half the total liquid water content (W), as follows:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{6}\times {\mathrm{10}}^{\mathrm{3}}}}& \underset{\mathrm{0}}{\overset{D_{\mathrm{0}}}{\int }}{D}^{\mathrm{3}}N\left(D\right)dD={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}}{\mathrm{6}\times {\mathrm{10}}^{\mathrm{3}}}},\\ & \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{D}^{\mathrm{3}}N\left(D\right)dD={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(W\right),\end{array}\end{array}$$

is also computed and included in the analysis.

Considering that a high-resolution dual-polarization X-band radar network is being deployed in Beijing for urban hydrometeorological applications, a series of polarimetric radar variables are simulated at X-band frequency based on the DSD measurements using the T-matrix method (Waterman, 1965; Leinonen, 2014), including horizontal reflectivity Z_H (mm⁶ m⁻³), differential reflectivity Z_dr (dB), and specific differential phase K_dp (^∘ km⁻¹). The drop-shaped model used in the simulation is the one proposed by Thurai et al. (2007). The temperature data are obtained from an automatic weather station collocated with the THUD disdrometer. In addition, various DSD-based radar QPE relations are derived and their parameterization errors are investigated for future development of the Beijing urban radar rainfall system.

2.3 Quality control

To minimize the measurement errors and improve data reliability, several quality control procedures have been applied to the 1 min DSD data. First, because of the low signal-to-noise ratios, the lowest two diameter bins are not used. That is, the raindrops less than 0.312 mm are eliminated in the analysis. Second, the 1 min sample data with total raindrop number smaller than 10 or the derived rain rate less than 0.1 mm h⁻¹ are considered noise and are removed (Sreekanth et al., 2017). Then, if the continuous data satisfying the above conditions last less than 5 min, they will be ignored to avoid spurious and erratic measurements (Jash et al., 2019). In addition, to focus on rainfall, all the data contaminated by hail are removed, and raindrops at a diameter of larger than 8 mm are eliminated (Bringi and Chandrasekar, 2001) since the biggest raindrops ever reported globally in the literature are around 8 mm (Baumgardner and Colpitt, 1995; Beard et al., 1986). Also, thresholds on the simulated radar parameters (i.e., $Z_{\mathrm{h}}=\mathrm{10}{\log }_{\mathrm{10}}{Z}_{\mathrm{H}}<\mathrm{55}$ dBZ, Z_dr>0 dB, and K_dp>0 ^∘ km⁻¹) are implemented to further guarantee the creditability of DSD data.

Figure 2 Histograms of different DSD parameters for all selected rainfall: (a) mass-weighted mean diameter, D_m (mm); (b) median volume diameter, D₀ (mm); (c) maximum diameter, D_max (mm); (d) generalized intercept parameter, log ₁₀N_w (N_w in m⁻³ mm⁻¹); (e) total number concentration, log ₁₀N_t (N_t in m⁻³); (f) rain rate, log ₁₀R (R in mm h⁻¹); (g) mass spectrum standard deviation, log ₁₀σ_m (σ_m in mm); (h) total number of raindrops, log ₁₀T_d; (i) liquid water content, log ₁₀W (W in g m⁻³).

3 DSD parameter characteristics

3.1 Distribution of DSD parameters

A total number of 43 618 1 min DSD spectra have been selected after data quality control, covering the wet seasons (May to October) from 2014 to 2018 except for May 2014 (no observation then). In this study, the raindrops below 1 mm are considered small drops; 1–3 mm are midsize drops; and large drops if larger than 3 mm (Krishna et al., 2016; Seela et al., 2017, 2018; Tokay et al., 2014). The distribution and statistics of the DSD parameters are shown in Fig. 2 and Table 1. D₀ and D_m have similar distributions, although D₀ has a larger range with a larger maximum and a smaller minimum value. It is more concentrated to small values, showing a smaller mean and median diameter with higher standard deviation, skewness and kurtosis values. The relationship $\mathrm{\Lambda }{D}_{\mathrm{m}}+\mathrm{3.67}=\mathrm{\Lambda }{D}_{\mathrm{0}}+\mathrm{4}$ (Ulbrich, 1983) may be explained for such a phenomenon when Λ>0. The distribution of D_max shows that during most of the rain events, the biggest drops are the middle class size, indicating that most of the rainfall is potentially made up of small and moderate raindrops. The statistical characteristics of log ₁₀N_w show almost equal median (3.596) and mean values (3.595), as well as a very small skewness value (0.040), indicating that log ₁₀N_w follows a symmetry distribution. The mean, median and skewness values of log ₁₀N_t, log ₁₀T_d, and log ₁₀σ_m also exhibit symmetry distributions. Moreover, the kurtosis of these three parameters is close to 3, which indicates that N_t, T_d, and σ_m obey the lognormal distribution. Since a threshold of 0.1 mm h⁻¹ is applied to the rain rate field (i.e., log ₁₀R is truncated by −1), the R meets a positive skew distribution. Because of this, log ₁₀W also has a positive skew distribution. It is worth noting that DSD samples with a rain rate about 0.8–1 mm h⁻¹ have the highest frequency and samples with a rain rate less than 1 mm h⁻¹ account for more than half of the total rain.

Table 1 Statistics of DSD parameters for all observations: D_m, D₀, D_max, log ₁₀N_w, log ₁₀N_t, log ₁₀R, log ₁₀σ_m, log ₁₀T_d and log ₁₀W.

3.2 DSD properties for different rain types

Previous studies in different climate regions have shown that DSD may substantially differ in the two general precipitation types (i.e., convective and stratiform), which has a great impact on the parameterization in both NWP models and remote sensing observations. In this study, rainfall events are separated into stratiform and convective cases using a method combining Bringi et al. (2003) and Chen et al. (2013). In particular, if the standard derivation of rain rate for a consequent 10 min is greater than 1.5 mm h⁻¹ and the rain rate is greater than 5 mm h⁻¹, it is classified as convective rain; otherwise, it is classified as stratiform rain.

Figure 3 shows the histograms of D_m and log ₁₀N_w for all the rainfall events and for the convective and stratiform subsets. The three key statistics are also indicted in Fig. 3, including mean, standard deviation (SD), and skewness. For the total dataset (Fig. 3a), the D_m histogram is highly positively skewed, while the skewness of log ₁₀N_w is near to zero, suggesting that the distribution of log ₁₀N_w is more symmetrical. The standard deviations of D_m and log ₁₀N_w are large (0.46 mm for D_m and 0.62 for log ₁₀N_w), indicating a high variability of both D_m and log ₁₀N_w. The mean values of D_m and log ₁₀N_w are 1.15 mm and 3.60, respectively. It should be noted that both mean values are slightly smaller compared with those obtained in the Beijing area during the summer time of 2015 (from 30 July to 30 September) and 2016 (from 9 June to 26 September) (G. Wen et al., 2017), which means that the DSD during summer time may be more concentrated than the whole rainy seasons.

Figure 3 Histograms of D_m and log ₁₀N_w for (a) all the rainfall events, (b) convective events, and (c) stratiform events. Mean values, standard deviation (SD), and skewness (SK) are also shown in the respective panels.

Considering different rain types, it can be found that the D_m for both types are positively skewed, while the skewness of log ₁₀N_w for convective is negative. The spread of log ₁₀N_w for convective rain is narrower compared to that of stratiform rain, and the skewness of log ₁₀N_w is larger than that of stratiform rain (−0.98 versus 0.10). The spreads and skewness of D_m for these two rainfall types perform oppositely (see Fig. 3b and c). In addition, histograms of D_m and log ₁₀N_w during convective rain tend to shift toward the large values relative to stratiform rain, indicating that convective events have higher D_m and log ₁₀N_w values than stratiform cases (1.91 mm and 3.66 for convective versus 1.08 mm and 3.59 for stratiform, respectively).

As Fig. 4 shows, in both convective and stratiform rains, with the increase in rain rate, the D_m increases (the positive exponents of the fitted power-law relationships), but the distributions of D_m become narrower. Note that at a higher rain rate, the D_m values tend to be stable, indicating that the DSD may have come to an equilibrium state where the coalescence and breakup of raindrops are in near balance (Hu and Srivastava, 1995). It can be seen in Fig. 4a that the D_m values reach a stable value around 2–2.5 mm when the rain rate R>60 mm h⁻¹, which means the increase in rain rate is mainly caused by an increase in concentration (Bringi and Chandrasekar, 2001). With respect to the D_m–R relationship, the coefficient and exponent values of convective rain are slightly higher than stratiform, suggesting a larger D_m of convective rain than stratiform rain for a given rain rate, which is different from the findings in eastern China (Wen et al., 2016) or southern China (A. S. Zhang et al., 2019).

Figure 4 Scatter density plot for D_m (mm) versus R (mm h⁻¹) for (a) convective events and (b) stratiform events. The fitted power-law relationships are also provided in each panel adopting a least-squares method.

Figure 5 shows the distribution of log ₁₀N_w versus D_m derived from the DSD data for all the rainfall events, as well as two different rain types. The two black rectangles correspond to the maritime and continental convective clusters and the dashed line corresponds to the stratiform (hereafter called “stratiform line”) case described by Bringi et al. (2003). For all the events combined, the distribution has a wide scale, but most points concentrate in the area of high log ₁₀N_w with low D_m. For convective and stratiform events, the distributions are concentrated in different areas (stratiform: 3.3–4.0 for log ₁₀N_w, 0.8–1.2 mm for D_m; convective: 3.7–4.2 for log ₁₀N_w, 1.4–2.0 mm for D_m). There is a rather clear boundary between the two rainfall types, although there are some overlaps. For convective rain, there are more points in the “Continental cluster” than the “Maritime cluster”, but most points are neither in the “Continental cluster” nor in the “Maritime cluster” and have a tendency to approach the stratiform rain. This indicates that the wet season convective rain in Beijing is neither maritime or continental as described by Bringi et al. (2003), which is likely due to the certain distance between Beijing and the nearest ocean (about 160 km). For stratiform rainfall, the points are more concentrated, even with a wide range of log ₁₀N_w versus D_m. More than 85 % of the stratiform points appear on the left side of the “stratiform line”. The average point of log ₁₀N_w−D_m for all the rainfall events combined (magenta hollow star) also appears on the left side of the “stratiform line” due to the highest population of stratiform in the summer monsoon season (see also Table 2). These indicate the lower diameter and higher concentration characteristics of rainfall in the Beijing area. The relationship of log ₁₀N_w−D₀ (see Fig. S1 in the Supplement) shows that the line to classify rain types based on log ₁₀N_w−D₀ (Thurai et al., 2016) would misclassify more convective rain as stratiform rain. This is probably due to the complex terrain in Beijing (Fig. 1a), where the high mountain to the west may have a substantial impact on the rain evolving from the western mainland.

Figure 5 Scatter density plot of log ₁₀N_w versus D_m: (a) the total rainfall events; (c) stratiform events; (d) convective events. (b) is the scatterplot of log ₁₀N_w versus D_m for the convective (red circle dots) and stratiform (blue square dots) cases. The two black rectangles in each subplot correspond to the maritime and continental convective clusters, and the black dashed line is the log ₁₀N_w–D_m relationship for stratiform rain reported by Bringi et al. (2003). The cross, hollow triangles, circle, squares, diamonds, and hearts in (b) represent the averaged values obtained in previous studies by Chen et al. (2013), Wen et al. (2016), G. Wen et al. (2017), Tang et al. (2014), and A. S. Zhang et al. (2019) for different parts of China. The colors of these symbols represent different events: magenta for total rainfall events; green for convective events; yellow for stratiform events; and black for the shallow events, a third type of precipitation besides convective and stratiform suggested by a few researchers, based on data from vertically pointing radar observations (Fabry and Zawadzki, 1995; Cha et al., 2009) in the study by Wen et al. (2016).

Table 2 Statistical properties of DSD parameters for convective and stratiform rain.

The comparison of DSDs in different parts of China shows interesting results. Even in the same region, the DSDs measured by different instruments have notable differences, such as the differences in Beijing between results from G. Wen et al. (2017) (2DVD, circle) and Tang et al. (2014) (Parsivel, square). In order to reduce the errors caused by different measurement instruments, only DSDs measured by Parsivel disdrometers are analyzed in this study. It is concluded that the eastern part of China has the lowest mean value of log ₁₀N_w (3.42) and highest mean value of D_m (1.66), while southern China has the highest mean value of log ₁₀N_w (3.86) with a middle value of D_m (1.46), and the northern part of China has the middle value of log ₁₀N_w (3.60) with a lowest value of D_m (1.15). This highlights that the DSD characteristics are highly dependent on the specific geographical locations and associated climate regimes. The results of Beijing from this study and Tang et al. (2014) show great differences in convective rain and lesser differences in stratiform rain, which is attributed to different convective systems during different years.

The DSD spectra and R(D) distributions of two rain types are shown in Fig. 6. Substantial differences are observed between these two rainfall types in both DSD spectra and R(D) distributions. The peaks of DSD spectra for both rainfall types are at the same diameter bin around D∼0.5 mm, while the spectrum for convective is higher than that of stratiform. The peak of R(D) distribution for stratiform rain is at the diameter around 0.9 mm and 1.9 mm for convective rain, which is much larger than where the DSD spectra peaks occur due to the D³ dependency of R(D). In addition, the distribution of R(D) for convective rain is much lower and broader. The differences in DSD spectra and R(D) distributions indicate that the convective rainfall has a higher concentration of moderate- to large-sized drops, and the large drops contribute more to convective rainfall compared to stratiform rainfall.

Figure 6 Composite raindrop spectra (a) and normalized R(D) distributions (b) for different rain types.

Figure 7 Same as Fig. 6 but for different rain rate classes.

3.3 DSD characteristics in different rain rate classes

To further understand the characteristics of DSD at different rainfall intensities in the Beijing area, the DSD measurements are divided into eight classes according to the associated rain rate (R): C1, $\mathrm{0.1}\le R<\mathrm{0.5}$; C2, $\mathrm{0.5}\le R<\mathrm{1}$; C3, $\mathrm{1}\le R<\mathrm{2}$; C4, $\mathrm{2}\le R<\mathrm{5}$; C5, $\mathrm{5}\le R<\mathrm{10}$; C6, $\mathrm{10}\le R<\mathrm{25}$; C7, $\mathrm{25}\le R<\mathrm{50}$; C8, R≥50 mm h⁻¹. Such classification is based on the fact of high frequency of low rain rates in the Beijing area as well as several previous studies, including Das and Maitra (2016), Harikumar et al. (2010), Krishna et al. (2016), Sarkar et al. (2015), and Tokay and Short (1996). The DSD sample numbers and rain rate statistics for each category are summarized in Table 3. For each rain rate class, the composite DSD spectrum is shown in Fig. 7a. Note that for almost all raindrop size bins, the concentration of a higher rain rate class is higher than that of a lower rain rate class. Furthermore, the breadth of the DSD shape increases and the tail of DSD shifts gradually to the larger diameter as the rainfall intensity increases, which is similar to previous findings in Taiwan (Seela et al., 2017), southern India (Jash et al., 2019), Palau (Krishna et al., 2016), and the UK (Islam et al., 2012). All the DSD spectra only have one peak, which differs from Krishna et al. (2016), where the spectrum becomes bimodal when the rain rate R>8 mm h⁻¹. In addition, the peaks of all DSD spectra are at a diameter around D∼0.5 mm, which is different from Jash et al. (2019) for India, where the peak position shifts towards larger diameters as the rain rate increases.

Table 3 Number and DSD retrieved rain rate statistics of each rain rate class.

The mean normalized R(D) of each rain rate class is shown in Fig. 7b, illustrating the contribution of each diameter bin to the total rainwater. The normalized rain rate distributions are unimodal and the peaks are around D∼0.9–2.5 mm. The peak position shifts to a larger diameter and the distribution becomes lower and broader as rain rate increases. These results are similar to those in Jash et al. (2019) for India but different from those in Peters et al. (2002) for Germany, where the R(D) distribution has a secondary peak at lower rain rate intensity (R<1 mm h⁻¹). This analysis implies that raindrops of diameter 0.9–2.5 mm (i.e., moderate size) contribute most towards accumulated rainwater during the rainy season in the Beijing area, and the size of drops contributing the most rainfall increases as the rainfall intensity increases.

Variations of the normalized intercept parameter (log ₁₀N_w) and mass-weighted mean diameter (D_m) in each rain rate class are provided in Fig. 8 with a box-whisker plot. It can be seen that D_m values increase with the increase in rainfall intensity, while the increasing trend of log ₁₀N_w is not as clear. This could be due to the imbalance between the decrease in small drop concentration and the increase in midsize and large drop concentration at a higher rain rate (R>10 mm h⁻¹, from C6 to C8). The means and standard deviations of D_m, log ₁₀N_w, N_t, W, μ, and Λ for each rain rate class are provided in Table 4, which clearly shows that with the increase in rainfall intensity, the mean values of total number concentration (N_t) and liquid water content (W) increase, while the mean values of the shape parameter (μ) and slope parameter (Λ) show a decreasing trend, resulting in a wider breadth and lower peak of DSD at high rain rates.

Figure 8 Variation of the normalized intercept parameter log ₁₀N_w (a) and the mass-weighted mean diameter D_m (b) for different rain rate classes. The white central line of the box indicates the median, the black central line in the box indicates the mean values, and the bottom and top lines of the box indicate the 25th and 75th percentiles, respectively. The bottom and top lines of the vertical lines out of the box indicate the 5th and 95th percentiles, respectively.

Table 4 Number and DSD retrieved rain rate statistics of each rain rate class.

3.4 Diurnal variations of DSD characteristics

Since the 1980s, Beijing has been experiencing rapid urbanization, causing a lot of problems, among which UHI is one of the most well-known phenomena (Yang et al., 2013b). Some studies showed that extreme precipitation events are more likely to occur during the period when the UHI intensity is high, usually from late afternoon to early morning in Beijing local time (LST) (Li et al., 2008; Song et al., 2014; Yang et al., 2013a, 2017; Y. Y. Zhang et al., 2019). In order to explore the DSD variations during the day, the diurnal periods are divided into four parts based on the UHI variation described in Yang et al. (2013b): strong UHI stage (S UHI, 21:00–06:00 LST), weak UHI stage (W UHI, 11:00–16:00 LST), UHI down stage characterized by a fast decline of UHI intensity (UHI D, 06:00–11:00 LST) and UHI up stage characterized by an abrupt rise of UHI intensity (UHI U, 16:00–21:00 LST). The rain rate and DSD characteristics corresponding to these four stages are shown in Table 5. The DSD spectra and R(D) distributions are shown in Fig. 9.

Figure 9 Same as Fig. 6 but for different diurnal periods based on UHI intensity.

Table 5 Mean and standard deviation (SD) values of R, D_m, log ₁₀N_w, N_t, W, μ, and Λ for different diurnal periods based on UHI intensity.

The DSD spectra of different diurnal periods are quite similar to those of different rain rate classes, showing a unimodal shape and peak position at the diameter D∼0.5 mm. It is notable that the DSD spectra are almost the same at small drop size bins (D<1 mm) and have the same width. As the diameter becomes larger, variations in the DSD spectra start showing up. The DSD spectra of the S UHI stage and UHI U stage show a similar and higher concentration, whereas the DSD spectra of the W UHI stage and UHI D stage have a similar but lower concentration, indicating that during the UHI U stage and S UHI stage, high-intensity rainfall is more likely to occur. This is in line with the study in Yang et al. (2017), which showed that the short-term high-intensity rainfall was more likely to happen at the UHI U stage and end at the late S UHI stage. The frequency and variation of the rain rate for different UHI stages (see Fig. S2) can also indicate this point.

The R(D) distributions for different diurnal periods in Fig. 9b show little difference between UHI U stage and S UHI stage, and the distributions at these two stages are lower and broader than the other two stages. At the W UHI stage, the R(D) distribution is the highest and the peak is at diameter around D∼0.9 mm, and the UHI D stage almost has the same peak around D∼0.9–1 mm, while the peaks during other two stages are at the diameter around D∼1 mm. That is, the drop size at the W UHI stage which contributes most to the accumulated rainwater is smaller than those at the UHI U stage or S UHI stage. The box-whisker plots of variation of D_m and log ₁₀N_w for each diurnal period show the same results (see Fig. 10). The W UHI stage has the highest mean concentration and the lowest mean D_m value, while the UHI U stage has the largest mean D_m value and the S UHI stage has the lowest mean concentration.

Figure 10 Same as Fig. 8 but for different diurnal periods based on UHI intensity.

3.5 DSD characteristics in different months

To obtain a better understanding of the seasonal variations of DSD characteristics in the Beijing urban area, rain data collected in different months are analyzed. The rain rate and DSD characteristics for different months are shown in Table 6. Figure 11 illustrates the corresponding DSD spectra and R(D) distributions.

Figure 11 Same as Fig. 6 but for different months.

Table 6 Mean and standard deviation (SD) values of R, D_m, log ₁₀N_w, N_t, W, μ, and Λ for each month.

As shown in Fig. 11, all the DSD spectra have peaks at diameter D∼0.5 mm, which are consistent with other classifications in this study. The DSD in May has a relatively higher concentration and a relatively lower concentration in July. At small drop size bins (D<1 mm), the spectra for May and September are similar, while the spectra for the other 4 months are similar. As the diameter increases, the differences between these spectra become larger, and the DSD spectrum for July has the highest concentration and October the lowest concentration. The rainfall with higher concentration and large drops is more likely to happen in July, leading to a high rain rate intensity (see also Fig. S3).

It is also noted that the R(D) distributions for each month are different from each other. The distributions of May, October, and September have a peak at diameter around D∼0.9 mm, while the distributions of June and August have a peak at diameter around D∼1 mm. The R(D) distribution of July has two peaks at diameter around D∼1 and D∼1.5 mm. In addition, the R(D) distribution of July is the widest and lowest, suggesting that a wide range of moderate drops contribute mostly to the rain in July. The D_m and log ₁₀N_w in Fig. 12 show an interesting annual cycle: the D_m (log ₁₀N_w) first goes up and (down) then goes down (up), while in July D_m (log ₁₀N_w) reaches the highest (lowest) value.

Figure 12 Same as Fig. 8 but for different months.

4 Implications for radar rainfall estimation

4.1 Single polarized radar applications

The power-law relationship between radar reflectivity (in mm⁶ m⁻³) and rain rate (in mm h⁻¹) (Z=aR^b) is the most widely used algorithm for single polarized radar QPE (including the current operational radars in Beijing). However, the coefficient a and exponent b greatly rely on the DSD variability (Bringi et al., 2003; Rosenfeld and Ulbrich, 2003; Uijlenhoet, 2001). The default Z–R relationship applied for the operational Weather Surveillance Radar – 1988 Doppler (WSR-88D) systems in the United States is Z=300R^1.4 (Fulton et al., 1998), whereas Z=200R^1.6 is commonly used in the continental area for stratiform rain (Marshall and Palmer, 1948, hereafter referred to as the MP-Stratiform relationship). The more appropriate and localized a and b are expected to improve regional radar rainfall estimation. In the following, the localized Z–R relationships for different rain types are derived by the nonlinear least-squares method, aiming to provide references for operational S-band radar rainfall applications in Beijing.

Figure 13 shows the scatter density plot of rain rate versus horizontal reflectivity, as well as the fitted power-law relations for different rain types. Figure 13 shows that most of the samples are at low values where both Z_H and R are small, which also suggests that the DSD may be under size-controlled conditions (Steiner et al., 2004). Meanwhile, the relationship for total rainfall (Z=238R^1.57) underestimates the rain rate at low values compared with the stratiform relationship (Z=171R^2.15), due to the inconsistent rain rate–reflectivity structures of two rain types.

Figure 13 Scatter density plot of R (mm h⁻¹) versus Z_H (mm⁶ m⁻³) for all rain events. The black, red, and blue curves, respectively, stand for the fitted power-law relations for total rain, convective rain, and stratiform rain. The purple and green dashed lines denote the default NEXRAD Z–R relation (Fulton et al., 1998) and a commonly used continental stratiform rain relation (Marshall and Palmer, 1948), respectively.

The default NEXRAD algorithm and MP-Stratiform relationship for continental stratiform rain are also indicated in Fig. 13 for comparison. At low reflectivity values (Z_H<23 dBZ), the curve of the MP-Stratiform relationship is below the local stratiform relation, but at higher values, it reverses. As the mean reflectivity of stratiform rain (21 dBZ) is less than 23 dBZ (see Table 2), the MP-Stratiform relationship may introduce underestimation of rainfall. The default NEXRAD relationship behaves similarly: underestimation at lower reflectivity values and overestimation at higher reflectivity values. Considering the mean reflectivity value of convective rain, the default NEXRAD relationship may cause overestimation of rainfall. In other words, the default relationship Z=300R^1.4 should be used with caution for local applications in Beijing.

4.2 High-frequency (X-band) polarimetric radar applications

A high-resolution dual-polarization X-band radar network is being deployed for urban hydrometeorological applications in the Beijing area. To support the radar deployment and facilitate the rainfall applications, the polarimetric parameters, including differential reflectivity Z_dr (dB) and specific differential propagation-phase shift K_dp (^∘ km⁻¹), are computed from the DSD measurements. Therein, the T-matrix method (Waterman, 1965) is adopted and the computations are made at X-band frequency. In addition, the polarimetric rainfall relations are derived based on the nonlinear least-squares method, including R(K_dp), R(K_dp, Z_DR), and R(Z_H, Z_DR). Here $Z_{\mathrm{DR}}={\mathrm{10}}^{Z_{\mathrm{dr}}/\mathrm{10}}$ is the differential reflectivity in the linear scale.

Figure 14 Scatter density plots of rainfall rates estimated from radar rainfall relations versus rain rates calculated directly from DSD: (a) R(Z_H), (b) R(K_dp), (c) (K_dp, Z_DR), and (d) R(Z_H, Z_DR). The black diagonal line in each panel represents the 1–1 relationship.

The derived X-band radar rainfall relations are as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }R\left(Z_{\mathrm{H}}\right)=\mathrm{0.0576}{Z}_{\mathrm{H}}^{\mathrm{0.557}},\text{(15)}& {\displaystyle }& {\displaystyle }R\left(K_{\mathrm{dp}}\right)=\mathrm{15.421}{K}_{\mathrm{dp}}^{\mathrm{0.817}},\text{(16)}& {\displaystyle }& {\displaystyle }R\left(K_{\mathrm{dp}},Z_{\mathrm{DR}}\right)=\mathrm{26.778}{K}_{\mathrm{dp}}^{\mathrm{0.946}}{Z}_{\mathrm{DR}}^{-\mathrm{1.249}},\text{(17)}& {\displaystyle }& {\displaystyle }R\left(Z_{\mathrm{H}},Z_{\mathrm{DR}}\right)=\mathrm{5.886}\times {\mathrm{10}}^{-\mathrm{3}}{Z}_{\mathrm{H}}^{\mathrm{0.994}}{Z}_{\mathrm{DR}}^{-\mathrm{4.929}}.\end{array}$$

Note that there are differences in the Z–R relationships between the X-band and S-band due to Mie scattering at higher frequency. Previous studies showed that the parameterization errors associated with various radar rainfall relations are among the key factors affecting the derived rainfall performance (Bringi and Chandrasekar, 2001). Hence, the parameterization errors in the X-band radar rainfall algorithms are investigated and quantified in this study. Figure 14 illustrates the scatter density plots of rain rates derived from R(Z_H), R(K_dp), R(K_dp, Z_DR), and R(Z_H, Z_DR) versus the rain rates directly computed from DSD. To quantify the parameterization errors, the normalized mean absolute error (NMAE) of the estimated rainfall rate is calculated, which is defined as

$$\begin{array}{}\text{(18)}& \mathrm{NMAE}={\displaystyle \frac{\langle |R_{\mathrm{EP}}-R_{\mathrm{D}}|\rangle }{\langle R_{\mathrm{D}}\rangle }},\end{array}$$

where the angle brackets stand for sample average; R_EP and R_D denote the estimated rain rates derived from parameterized radar rainfall algorithms and DSD information, respectively. The NMAE_RR is calculated for different rainfall rate intervals from 0 to 100 mm h⁻¹. Figure 15 shows the parameterization error structure of R(Z_H), R(K_dp), R(K_dp, Z_DR), and R(Z_H, Z_DR) as a function of rainfall rate.

It can be seen from Figs. 14 and 15 that the algorithms based on dual-polarization radar parameters can provide better estimates than the Z–R relationship. In addition, the dual-parameter algorithms, namely R(K_dp, Z_DR) and R(Z_H, Z_DR), have even better performance than the single-parameter-based algorithm including R(K_dp). The NMAE has a decreasing trend as the rain rate increases from 1 to 60 mm h⁻¹. The fluctuation when rain rate is greater than 60 mm h⁻¹ may be due to the random errors caused by a few samples of large values. The parameterization errors of R(K_dp), R(K_dp, Z_DR), and R(Z_H, Z_DR) become stable when rain rate becomes higher than 10 mm h⁻¹. It is also noted that at low rain rate (less than 10 mm h⁻¹), the NMAE of R(Z_H, Z_DR) is the smallest, while at higher rain rate (higher than 10 mm h⁻¹) the NMAE of R(K_dp, Z_DR) becomes the smallest. This again highlights the importance of selecting appropriate rain rate relations for local radar applications.

Figure 15 Parameterization error structure of R(Z_H), R(K_dp), R(K_dp, Z_DR), and R(Z_H, Z_DR) as a function of rainfall rate: (a) for mean rain rate less than 10 mm h⁻¹; (b) for rain rate of the whole dataset.

5 Summary and conclusion

In this paper, 5-year (2014–2018) observations of DSD from a disdrometer deployed at Tsinghua University are analyzed to explore the microphysical characteristics of precipitation during rainy seasons (May–October) in the Beijing urban area. The main conclusions are as follows.

1. For all rain events, all the DSD parameters (D_m, D₀, D_max, log ₁₀N_w, log ₁₀N_t, log ₁₀R, log ₁₀σ_m, log ₁₀T_d and log ₁₀W) except σ_m have a positive skewness, indicating a high frequency of low values and a low frequency of high values in the Beijing urban area. More than half of the DSD measurements are characterized by a rainfall rate of less than 1 mm h⁻¹.

2. The mean values of log ₁₀N_w and D_m of convective rain are higher than that of stratiform rain, indicating a higher raindrop concentration and larger drop size during convective events. This is also in line with the raindrop spectra and normalized R(D) distribution. In addition, log ₁₀N_w of convective rain is negatively skewed, which is opposite to that of stratiform rain. For both rainfall types, the D_m values are higher but the distributions are narrower at higher rainfall intensities.

3. There is a clear boundary to distinguish between convective and stratiform rain from the scatterplot of log ₁₀N_w versus D_m. However, the convective rain in the Beijing area is neither continental nor maritime as described by Bringi et al. (2003), due to the particular location and complex topography. Moreover, the comparison with different parts of China shows that the DSD variability is closely related to geographic location, climate regimes and study periods.

4. Stratified by rain rate, the DSD spectra become higher and wider as the rain rate increases, but all have peaks at the similar diameter size D∼0.5 mm. The peaks of the normalized R(D) distribution shift to larger diameter size (still within the midsize range) and the distribution becomes lower and wider as the rain rate increases. Meanwhile, the D_m and log ₁₀N_w show an increasing trend and the slope parameter (μ) shows a decreasing trend as the rain rate increases.

5. During the periods of strong UHI and UHI up stages, the DSD spectra trend to have a higher concentration at large size drops and larger D_m values than other periods, indicating intense rainfall during these periods. The DSD has similar characteristics in July and August. In addition, the D_m and log ₁₀N_w values show a diurnal cycle and an annual cycle. All these findings indicate substantial temporal variabilities of DSD in Beijing.

6. The Z–R relationship derived from local DSD in Beijing is quite different from the operational NEXRAD algorithm (MP-Stratiform) which may overestimate (underestimate) rainfall at high (low) rain intensity. The error structures of different algorithms show that the polarimetric radar rainfall relations R(K_dp), R(K_dp, Z_DR), and R(Z_H, Z_DR) have greater potential than Z–R methods for urban QPE.

The DSD characteristics presented in this study not only provide further understanding of precipitation microphysical variabilities in Beijing, but also provide indications of future model development to improve local precipitation forecast. In addition, a high-resolution X-band dual-polarization radar network is being deployed in Beijing. This study is expected to provide references for future development of localized radar rainfall algorithms. Nevertheless, the DSD spectra also show the limitations of the Parisvel² disdrometer in measuring small raindrops. Future study should be carried out with multiple instruments including a two-dimensional video disdrometer just deployed in this area. We also want to note that combining additional observations such as the vertically pointing profiler radar data (White et al., 2003) can further enhance the classification results of different rainfall types, which should be considered in future studies. In addition, further investigation on the spatial variability of DSD induced by the complex micro-topography in urban areas should be conducted.