Vertical profiles of raindrop size distribution parameters of summer rainfall in the eastern Tibetan Plateau: retrieval method and characteristics

The eastern Tibetan Plateau has a high elevation, with a cold and dry atmospheric background. The features of the raindrop size distributions (DSD) in this region have notable differences from those in the plains. The general empirical relationships for retrieving parameters of precipitation from radar observations are not applicable in the eastern Tibetan Plateau. In this study, we developed a new method based on optimal estimation theory to retrieve the vertical profiles of N₀ and D_m from a Ka-band zenith-pointing Doppler radar. Validation by a field campaign during the summer of 2024 indicates that the mean bias in the log ₁₀(N₀) and D_m derived from the PARSIVEL² disdrometer and the retrieved values are 0.12 and −0.1 mm respectively, demonstrating the effectiveness of the retrieved DSD parameters in this region. Based on the retrieved vertical profiles of DSD parameters, some unique characteristics are found. The heavy precipitation (the maximum value in the reflectivity profile exceeding 30 dBZ) exhibits a higher particle number concentration above 2 km and larger raindrop size in the bottom of the rainfall on average. The mean values of D_m above 2 km are approximately 0.5 mm, for heavy precipitation, the value increase as the raindrops fall, reaching a peak at around 0.5 km. Precipitation that occurs after the nighttime cooling usually has higher particle concentrations and smaller particle sizes. Based on the above research, empirical relationships for the quantitative precipitation estimates (QPE) and attenuation correction using Ka-band radar in the eastern Tibetan Plateau are established.

1 Introduction

Modeling and forecasting of precipitation is one of the key challenges in the field of meteorology, with substantial implications for quantifying the dynamics of water and energy cycles within the Earth-atmosphere system and improving the accuracy of numerical weather and climate models (Maggioni and Massari, 2019; Foufoula-Georgiou et al., 2020). The parameters of raindrop size distributions (DSD) are critical for describing the microphysical processes during precipitation (Thomas et al., 2021). Exploring the dynamic evolution of the vertical profile of DSD parameters may enhance the understanding of precipitation processes (Tapiador et al., 2010). Also, DSD is associated with the backscattering signals of radar observations, thereby affecting the outcomes of quantitative precipitation estimates (QPE) derived from radar (Gatidis et al., 2022). Furthermore, DSD and their associated uncertainties are critical components of data assimilation systems and numerical weather models (Yang et al., 2019).

The eastern Tibetan Plateau serves as a pivotal area for the transport of moisture from the Indian Ocean to eastern Asia (Xu et al., 2019). The water cycle in this region affects precipitation patterns across China and the global climate (Hu et al., 2021). The average altitude of this region is above 3000 m, and the atmospheric density is low, with a sparse water vapor content (You et al., 2015). The unique atmospheric background in this region results in distinct features of DSD compared to plain areas (Li et al., 2022). Exploring the characteristics of the vertical distribution of DSD parameters in the eastern Tibetan Plateau is crucial for improving the accuracy of local microphysical parameterizations of precipitation and deepening the understanding of the water cycle in this region.

The high temporal and spatial resolution of weather radar can provide continuous profiles of the precipitation parameters. However, the precipitation intensity and the size of raindrops in the eastern Tibetan Plateau are relatively small compared to those in the middle and lower reaches of the Yangtze River (Pu et al., 2021), and the sensitivity of common weather radars to small particles is limited at their operating frequencies. Ka-band radar exhibits superior spatial resolution and sensitivity to small precipitation particles, and has been widely applied in the observation of precipitation microphysical parameter profiles (Clothiaux et al., 2007). Studies have indicated that the attenuation of the reflectivity factor (Z) of Ka-band radar caused by large raindrops cannot be ignored (Daisley and Aydin, 2002). To minimize the impact of the attenuation effect on the retrieval accuracy, the Doppler spectrum observed by Ka-band radar was used to retrieve the parameters of DSD (Tridon and Battaglia, 2015; Giangrande et al., 2012). However, the results may be affected by turbulence of air, and the retrieval may have multiple solutions (Williams et al., 2016). Observations of the multi-frequency radar can provide additional information on the parameters of DSD. The dual-frequency ratio and differential attenuation coefficients were used to constrain the retrieval of DSD parameters (Unal and Van Den Brule, 2024; Hogan et al., 2005). Yet, the empirical relationships between the precipitation parameters and the observations of radar were affected by the average features of DSD that were summarized from plain areas, and are therefore not applicable in the eastern Tibetan Plateau (Protat et al., 2019). Further research is needed to accurately retrieve the profiles of DSD parameters from Ka-band radar in the eastern Tibetan Plateau and to summarize the characteristics of the vertical distribution of DSD in this region.

Due to the limited information obtained from single-frequency radar, the methods for retrieval of DSD parameters are generally chosen to fix the μ parameter of the Gamma Size Distribution (GSD) (Kumar et al., 2011) or utilize λ-μ empirical relationships of GSD to constrain the retrieval results (Huang et al., 2021), where μ and λ are shape, and slope parameters of GSD, respectively. However, these assumptions inevitably introduce errors into the retrieval results. Although the GSD can provide a detailed structure of the DSD, previous studies have indicated that the Marshall–Palmer exponential (MP) drop size distribution is sufficient to represent the size distribution of particles for cloud liquid water and precipitation (Smith, 2003). Besides, the representation of rainfall in the weather and climate models usually adopts a single or two-parameter scheme for rainfall to reduce the complexity of calculation (Jin and Baik, 2023). This article presents a retrieval method developed for Ka-band zenith-pointing radar and the characteristics of retrieved vertical profiles of DSD parameters. The retrieval method is based on optimal estimation theory, corrects the precipitation attenuation layer by layer and retrieves the N₀ and D_m vertical profiles of MP distribution, and provides uncertainty estimates for the retrieval. This study primarily addresses the precipitation scenarios where the melting layer can be observed in radar LDR (Linear Depolarization Ratio) measurements. The retrieval and characteristics of DSD parameters under strong convective conditions are not included in this study. In Sect. 2, the instruments and retrieval method used in this study are introduced. Section 3 presents the validations of the forward calculation and retrieval method. The spatiotemporal characteristics of the vertical profiles of the DSD parameters and the fitted empirical relationships are given in Sect. 4.

2 Instrument and retrieval method

During the 2024 summer precipitation observation experiment in the eastern Tibetan Plateau, a comprehensive study of cloud-precipitation was conducted at the Hongyuan Meteorological Bureau weather site (32.8° N, 102.5° E, 3499 m above sea level). In this study, the vertical distribution of DSD parameters of summer rainfall in the eastern Tibetan Plateau is investigated using a Ka-band zenith-pointing Doppler radar and a ground-based disdrometer deployed at the Hongyuan site. The distance between the two instruments is about 70 m.

2.1 Ground-based disdrometer

The surface DSD are collected by the PARSIVEL² disdrometer. The measuring ranges of the disdrometer are 0.062–24.5 mm for the size of raindrops and 0.05–20.8 m s⁻¹ for velocity, respectively. The time resolution is 60 s. According to the number of raindrops in the discrete scale intervals observed by the disdrometer n and the terminal falling velocity V, the raindrop size distribution N(D) in ${\mathrm{mm}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{3}}$ is written as:

$$\begin{array}{}\text{(1)}& N\left(D_i\right)=\sum _{j=\mathrm{1}}^{\mathrm{32}}{\displaystyle \frac{n_{i,j}}{A_i\cdot t\cdot V_{i,j}\cdot \mathrm{\Delta }{D}_i}},\end{array}$$

Where n_i,j represents the number of raindrops in the i-th size bin and the j-th velocity bin, A_i denotes the sampling area, t is the integration time, and V_i,j indicates the falling velocity of raindrops in the i-th size bin and the j-th velocity bin. The particle diameter is expressed in millimeters (mm), while the velocity is measured in meters per second (m s⁻¹) in this study. The MP exponential drop size distribution (Smith, 2003) expressed as:

$$\begin{array}{}\text{(2)}& N\left(D\right)=N_{\mathrm{0}}\exp (-\mathrm{\Lambda }D)=N_{\mathrm{0}}\exp (-\mathrm{4}D/D_{\mathrm{m}}),\end{array}$$

In Eq. (2), D is the maximum dimension of the particle, N₀ is the concentration parameter, the mass-weighted mean diameter D_m, is determined as the ratio of the fourth to third moments of the DSD. N₀ is calculated from the third and fourth moments of DSD (Smith, 2003), expressed as:

$$\begin{array}{}\text{(3)}& N_{\mathrm{0}}=\mathrm{128}{M}_{\mathrm{3}}^{\mathrm{5}}/\left(\mathrm{3}{M}_{\mathrm{4}}^{\mathrm{4}}\right)\end{array}$$

The rain rate R is defined by the fall rate and DSD as:

$$\begin{array}{}\text{(4)}& R=\int n\left(D\right)(\mathit{\pi }{D}^{\mathrm{3}}/\mathrm{6})V\left(D\right)\mathrm{d}D\end{array}$$

For Ka-band radar, the reflectivity factor Z and the Doppler radial velocity V_r are calculated from Eqs. (5) and (6).

$$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{rl}Z& ={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K{|}^{\mathrm{2}}}}\sum _{i=\mathrm{1}}^{\mathrm{32}}N\left(D_i\right){\mathit{\sigma }}_{\mathrm{b}}\mathrm{\Delta }{D}_i\\ & ={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{4}}|K{|}^{\mathrm{2}}}}\sum _{i=\mathrm{1}}^{\mathrm{32}}N\left(D_i\right)Q_{\mathrm{b}}\left(D\right)D^{\mathrm{2}}\mathrm{\Delta }D\end{array}\text{(6)}& {\displaystyle }{V}_{\mathrm{r}}\left(D\right)={\displaystyle \frac{\underset{D_{\text{min}}}{\overset{D_{\text{max}}}{\int }}N\left(D\right)V\left(D\right){\mathit{\sigma }}_{\mathrm{b}}\mathrm{d}D}{\underset{D_{\text{min}}}{\overset{D_{\text{max}}}{\int }}N\left(D\right){\mathit{\sigma }}_{\mathrm{b}}\mathrm{d}D}}\end{array}$$

Where the efficiency factor of scattering (Q_b), backscattering section of radar (σ_b), and $|K{|}^{\mathrm{2}}$ are calculated by pyQuickBeam (Haynes et al., 2007). In Eq. (6), V(D) is the final falling velocity of raindrops considering the correction term of altitude h, which is calculated from Eq. (7) of Atlas et al. (Atlas et al., 1973):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\displaystyle }D& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{0.6}}}\times \ln {\displaystyle \frac{\mathrm{10.3}}{\mathrm{9.65}-V\left(D\right)/\mathit{\delta }\left(h\right)}}\\ {\displaystyle }\mathit{\delta }\left(h\right)& {\displaystyle }=\mathrm{1}+\mathrm{3.68}\times {\mathrm{10}}^{-\mathrm{5}}h+\mathrm{1.71}\times {\mathrm{10}}^{-\mathrm{9}}{h}^{\mathrm{2}}\end{array}\end{array}$$

Historical observations of the disdrometer at Hongyuna site from June to August of 2021 to 2023 are used to obtain prior information on the parameters of the MP size distribution of summer rainfall in this region. Based on the observed DSD, the corresponding Z of Ka-band radar is calculated using Eq. (5). The calculated Z values are categorized into four intervals: $[-\mathrm{10},\mathrm{10})$, [10,20), [20,30), and [30,45). The log ₁₀(N₀), D_m and Rain Water Content (RWC) in each Z interval are counted separately and illustrated in Fig. 1. In the box plot used in this paper, the box indicates the Interquartile Range (IQR, 25th–75th percentiles), the horizontal line inside the box denotes the median, and the whiskers extend to the most extreme values within 1.5×IQR. Circles denote outliers beyond the whiskers. RWC is calculated by, $\text{RWC}=\frac{\mathit{\pi }}{\mathrm{6}}\times {\mathrm{10}}^{-\mathrm{4}}{\int }_{D_{\text{min}}}^{D_{\text{max}}}N\left(D\right)D^{\mathrm{3}}\mathrm{d}D$. Approximately 113 000 samples of DSD are collected. The results show that the distribution of the log ₁₀(N₀) ranges from 1.0 to 5.6. The majority of D_m is less than 1.7 mm. The statistical results indicate that as the Z of precipitation increases, the proportion of larger raindrops of DSD rises, and the concentration of raindrops decreases. Since the RWC is related to the 3-th moment of DSD, it increases as the rain intensity increase and exceeds on average about 0.2 g m⁻³ at strong echoes (30–45 dBZ).

Figure 1 Box plots of the parameters of MP size distribution in different Z intervals. (a) log ₁₀(N₀). (b) D_m. (c) RWC.

2.2 Ka-band zenith-pointing observation radar

The Ka-band zenith-pointing radar operates at the frequency of 35 GHz. The observation uncertainties in Z, V_r, and LDR are 0.5 dBZ, 0.5 m s⁻¹, and 0.5 dB, respectively. The detection range of radar is 0.12 to 20 km, with a vertical resolution of 30 m and a beam width less than 0.4°. The integration time is 60 s. To mitigate the impact of near-field uncertainties of radar observations, the seventh range gate (180–210 m) of the radar is chosen as the first effective observation at the surface, denoted as h₀. The “0” of h₀ refers to the index of the range gate of the effective observations of radar. When the Z measured by radar at h₀, $Z_{h_{\mathrm{0}}}$, exceeds −10 dBZ, it is identified as a precipitation scenario. The range gate of the precipitation height (h_top) is determined from the bottom of the melting layer assessed by the LDR observations. The method (Thurai et al., 2016) for determining the htop is as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\displaystyle }{\text{LDR}}_{h_{\text{top}}}& {\displaystyle }>max\left(\right[{\text{LDR}}_{h_{\mathrm{0}}},\mathrm{\dots },{\text{LDR}}_{h_{(\text{top}-\mathrm{1})}}\left]\right)\\ {\displaystyle }{\text{LDR}}_{h_{\text{top}}}& {\displaystyle }-{\text{LDR}}_{h_{(\text{top}-\mathrm{2})}}\ge \mathrm{10}\mathrm{dB}\end{array}\end{array}$$

After determining the height of the precipitation layer, the parameters of the DSD are retrieved from bottom to top. For each range gate of the radar observation, the optimal estimation method is used to retrieve log ₁₀(N₀) and D_m from Z and V_r.

The optimal estimation is a widely used physical retrieval method that iterates the maximum probability state x_op of the parameters to be retrieved by combining observations, prior information, and associated uncertainties (Löhnert and Maier, 2012). The optimal estimation algorithm employed in this study is based on the publicly available code by Maahn et al. (Maahn et al., 2020). The observational vector y_obs of the optimal estimation algorithm consist of the Z and V_r observed by the Ka-band zenith-pointing radar, and the state vector x comprises log ₁₀(N₀) and D_m. The first guess of x_a, and a priori error covariance matrix S_a are set based on the prior information of the DSD in Hongyuan. From the precipitation cases described in Sect. 2.1, a subset of 6680 case is selected to construct the prior dataset. The dataset is obtained by stratified random sampling based on the calculated Z. This approach reduces the over representation of frequent weak precipitation samples in the prior and yields a more uniformly distributed and representative priori dataset (Boschetti et al., 2016). Then, the values of D_m are normalized since the state vector are assumed to be approximately Gaussian in optimal estimation, and statistical analysis of the dataset shows that the distribution of D_m departs from normality. After the transformation, the prior database for the retrieval is established. The mean value of log ₁₀(N₀) and D_m in the dataset is set as the first guess x_a for the retrieval. The frequency distribution of log ₁₀(N₀) and normalised D_m are plotted in Fig. 2a. S_a calculated from the database are plotted in Fig. 2b. In the study, the S_a is a constant matrix, and includes correlations between N₀ and D_m derived from the prior database.

Figure 2 The frequency distribution of the prior database and the heatmap of S_a used in the retrieval. (a) Frequency distribution of N₀ and normalized D_m, and (b) heatmap of S_a.

The error covariance matrix S_y is set according to the observation uncertainty of Z and V_r. The measurement of Z and V_r is assumed to be uncorrelated, and the value on the diagonal of the matrix is squire of the observation uncertainty. The x_a is taken as the input for Eqs. (5) and (6) to calculate the simulated observation F(x_a). The S_y and S_a are then used to weight the difference of the F(x_a) and y_obs for updating the x_a to x_j+1. The subscript is the iteration index. X_j+1 is calculated by:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\displaystyle }{\mathit{x}}_{j+\mathrm{1}}& {\displaystyle }={\mathit{x}}_{\mathrm{a}}+({\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}+{\mathbf{K}}_j^{\mathrm{T}}{\mathbf{S}}_y^{-\mathrm{1}}{\mathbf{K}}_j{)}^{-\mathrm{1}}\\ {\displaystyle }& {\displaystyle }{\mathbf{K}}_j^{\mathrm{T}}{\mathbf{S}}_y^{-\mathrm{1}}[\mathit{y}-F({\mathit{x}}_j)+{\mathbf{K}}_j({\mathit{x}}_j-{\mathit{x}}_{\mathrm{a}}\left)\right]\end{array}\end{array}$$

Where K_j represents the Jacobian matrix computed at the j-th iteration. The elements of K_j numerically computed by finite differences around the state vector at each iteration of the retrieval, the perturbation is chosen as 1 % of the state vector. The typical value of K_j indicate that the Z is sensitive to both log ₁₀(N₀) and D_m, with $\partial Z/\partial {\log }_{\mathrm{10}}\left(N_{\mathrm{0}}\right)$ and $\partial Z/\partial D_{\mathrm{m}}$ on the order of 10 and 20, respectively. And V_r is almost insensitive to log ₁₀(N₀) and only affected by D_m, with $\partial V_{\mathrm{r}}/\partial {\log }_{\mathrm{10}}\left(N_{\mathrm{0}}\right)\approx \mathrm{0}$ and $\partial V_{\mathrm{r}}/\partial D_{\mathrm{m}}\approx \mathrm{2.3}$. This pattern indicates that the mainly constrains the retrieval of N₀ and D_m, while V_r provides complementary sensitivity to D_m. The iteration ceases when the result of the left side of Eq. (10) falls below the predetermined threshold χ.

$$\begin{array}{}\text{(10)}& ({\mathit{x}}_{j+\mathrm{1}}-{\mathit{x}}_j{)}^{\mathrm{T}}{\mathbf{S}}_j^{-\mathrm{1}}({\mathit{x}}_{j+\mathrm{1}}-{\mathit{x}}_j)=\mathit{\chi }\end{array}$$

χ is typically set to the length of the parameters that are to be retrieved. S_j in Eq. (11) provides the uncertainty associated with the retrieved x_j:

$$\begin{array}{}\text{(11)}& {\mathbf{S}}_j=({\mathbf{S}}_j^{-\mathrm{1}}+{\mathbf{K}}_j^{\mathrm{T}}{\mathbf{S}}_y^{-\mathrm{1}}{\mathbf{K}}_j{)}^{-\mathrm{1}}\end{array}$$

If the convergence of retrieval is achieved, the x_j and S_j are the optimal solution x_op and corresponding uncertainty S_op.

The retrieval process of the vertical profiles of the DSD parameters is described as follows. The temperature at $Z_{h_{\mathrm{0}}}$ is set to 287.15 K based on statistical results of temperature from the microwave radiometer at a height of 210 m in Hongyuan. Unless otherwise noted, heights throughout this paper refer to above ground level (AGL) height. The temperatures at different range gates of the radar are calculated using the temperature lapse rate. The temperature lapse rate is set to 6.5 ° km⁻¹ (Kattel et al., 2015). The $Z_{h_{\mathrm{0}}}$ is assumed to be unaffected by the attenuation effect. The $Z_{h_{\mathrm{0}}}$ and $V_{h_{\mathrm{0}}}$ are taken as the y_obs for the optimal estimation algorithm to retrieve the log ₁₀(N₀) and D_m of the DSD at the surface. Then, the N₀ and D_m at the surface are used to correct the Z of the next range gate $Z_{h_{(i+\mathrm{1})}}$:

$$\begin{array}{}\text{(12)}& Z_{h_{(i+\mathrm{1})}}=Z{a}_{h_{(i+\mathrm{1})}}\exp \left[-\mathrm{2}\underset{\mathrm{0}}{\overset{h_i}{\int }}k\left(h\right)\mathrm{d}h\right]\end{array}$$

where $Z_{h_{(i+\mathrm{1})}}$ and $Z{a}_{h_{(i+\mathrm{1})}}$ are the attenuated and actual radar reflectivity in mm⁶ m⁻³, k is the specific attenuation in dB km⁻¹, defined as:

$$\begin{array}{}\text{(13)}& k=\underset{D_{\text{min}}}{\overset{D_{\text{max}}}{\int }}{\mathit{\sigma }}_{\mathrm{e}}N\left(D\right)\mathrm{d}D=\underset{D_{\text{min}}}{\overset{D_{\text{max}}}{\int }}{\displaystyle \frac{\mathit{\pi }{D}^{\mathrm{2}}}{\mathrm{4}}}{Q}_{\mathrm{e}}N\left(D\right)\mathrm{d}D\end{array}$$

Q_e and σ_e is calculated from the drop size and temperature using the pyQuickBeam software. If the height of the following range gate, h_(i+1), is less than htop, the observational vector of the next range gate is updated to ${\mathit{y}}_{\text{obs}\mathit{\_}h(i+\mathrm{1})}=[Z{a}_{h_{(i+\mathrm{1})}},V_{h_{(i+\mathrm{1})}}]$. Next, the optimal estimation is used to retrieve the log ₁₀(N₀) and the D_m at h_(i+1). Following this procedure, the observations of Z are corrected layer by layer, and the N₀ and D_m are retrieved bottom to top. The flowchart of the retrieval method is presented in Fig. 3.

Figure 3 Flowchart for the retrieval method of vertical profiles of DSD parameters.

3 Validation

3.1 Validation of the forward calculation

The accuracy of forward calculations is crucial for the optimal estimation algorithm. The observations of the Ka-band zenith-pointing radar and the disdrometer during the field campaign are used to validate the forward calculations. The Z and V_r observed by the radar at h₀ serve as the truth values, while the N(D) observed from the PARSIVAL² disdrometer is used as the input to Eqs. (5) and (6) to calculate the simulated Z and V_r. The mean error (ME), Pearson correlation coefficient (Pearson's r), and Jensen-Shannon divergence (JSD) are employed to assess the bias between the true values and the simulated values. The ME is used to evaluate the average difference between two datasets. Pearson's r measures the linear relationship between the two datasets. The JSD quantifies the similarity in the probability distributions of two datasets (Zhang et al., 2023).

Figure 4 illustrates the true values of Z and V_r observed from the Ka-band zenith-pointing radar and simulated values from the disdrometer during the field campaign in July and August 2024. According to the assessment metrics, the true values and the simulated values of Z and V_r are in good agreement. The ME of Z and V_r is 3.48 dBZ and $-\mathrm{1.04}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The Pearson's r of Z and V_r are above 0.88. For precipitation scenarios with $Z_{h_{\mathrm{0}}}<\mathrm{25}\mathrm{dBZ}$, the Z calculated from the ground-based disdrometer observations generally exceeds the radar measurements. For radar observations with Z≥25 dBZ, the reflectivity factors derived from disdrometer measurements are typically lower than those observed by radar. The simulations and observations for V_r less than 3 m s⁻¹ are consistent, while for V_r exceeding 3 m s⁻¹, the simulations tend to underestimate the values of V_r. JSD values of the simulated and observed suggest some mismatch in Z and V_r distributions. This bias is possibly due to differences in the sampling volumes of the radar and disdrometer, since the disdrometer measures the raindrop at surface and the first effective range bin of radar is at 210 m. In addtion, the wind at the surface may influence the accuracy of the measurements of the disdrometer. Furthermore, as the raindrops fall, they may deviate from the spherical shape assumption in the radar forward model, leading to discrepancy between the forward simulation and observed values.

Figure 4 The true values of Z and V_r observed from the Ka-band zenith-pointing radar and simulated values during the field campaign in July and August 2024. (a) Z. (b) V_r.

3.2 Validation of the optimal estimation algorithm

To validate the effectiveness of the optimal estimation algorithm, a random stratified sampling of 1800 cases of DSD is selected from the historical observations of the ground-based disdrometers. The precipitation cases were stratified according to the distribution range of N₀ and D_m, and 900 cases were sampled for each of the two parameters to ensure that the validation data covered a relatively comprehensive precipitation scenario. Random noise within 0.5 dBZ and 0.5 m s⁻¹ is added to the values of Z and V_r, which are calculated from the selected DSD cases with Eqs. (5) and (6). Then the Z and V_r are taken as the input for the optimal estimation algorithm to retrieve the log ₁₀(N₀) and D_m of DSD. The log ₁₀(N₀) and D_m derived from the selected DSD with Eq. (3) are taken as the true values of the retrievals. The final value of S_x after retrieval convergence is indicated by error bar in the figure. The comparison of the retrieved values and the true values of log ₁₀(N₀) and D_m is presented in Fig. 5. The average bias of the retrieved log ₁₀(N₀) and D_m is −0.056 and 0.015 mm, with the average uncertainties of 0.41 and 0.15 mm. The correlation between the retrieved values and the true values of log ₁₀(N₀) and D_m both exceeds 0.9. The JSD values of 0.08 and 0.16 indicate a good consistency between the distributions of the retrieved values and the true values of log ₁₀(N₀) and D_m. The source of the retrieval errors may be attributed to the random noise added to the Z and V_r. Furthermore, the errors of the retrieval are primarily observed in scenarios where the D_m values exceeded 2.5 mm and the log ₁₀(N₀) value is less than 2.6. This may be attributed to the lack of prior information in these cases.

Figure 5 The retrieved values and the true values of log ₁₀(N₀) and D_m. (a) log ₁₀(N₀). (b) D_m.

Figure 6 The retrieved and observed log ₁₀(N₀) and D_m during the field campaign in July and August 2024. (a) log ₁₀(N₀). (b) D_m.

The optimal estimation algorithm is further validated with the observations from the Ka-band zenith-pointing radar and the ground-based disdrometer in Hongyuan during the field campaign in July and August 2024. The Z and V_r observed by radar at h₀ are used as inputs for the optimal estimation algorithm. The log ₁₀(N₀) and D_m derived from the disdrometer are used for the comparison. Approximately 5900 effective observations of the radar and the disdrometer of precipitation are used for the validation. In Fig. 6, the retrieved values from the optimal estimation algorithm are denoted with black points, while the uncertainties of the retrieval are indicated by error bars in gray. The parameters of DSD derived from the disdrometer observations are marked with blue points. The correlation coefficients between the retrieved and the derived log ₁₀(N₀) and D_m are 0.63 and 0.76. The ME of the retrieved and the derived log ₁₀(N₀) and D_m are 0.12 and −0.01 mm, and average uncertainties are 0.41 and 0.16 mm, respectively. Most of the log ₁₀(N₀) and D_m values derived from the disdrometer are within the uncertainty range of the retrievals. The JSD values indicate that the retrieved D_m values are consistent with the distribution of the observations of the ground-based disdrometer. The bias between the retrieved and observed parameters of DSD may be attributed to the observational uncertainties and the differences in the sampling volumes of the two instruments, as well as errors introduced when fitting the DSD parameters. More importantly, based on the validation of the radar forward calculation, there is a bias between the radar observation and simulation. Physical assumptions in the scattering calculations, such as the spherical-drop approximation, may also contribute to discrepancies between the retrieved values and the surface DSD parameters.

Figure 7 Profiles of the attenuation-corrected Z, the retrieved profiles of DSD parameters, and PLA of the precipitation on 15 August 2024. (a) Z. (b) log ₁₀(N₀). (c) D_m and (d) PIA.

4 Results

4.1 Case analysis

Figure 7 presents the attenuation-corrected profiles of Z for the precipitation on 15 August 2024, along with the retrieved vertical profiles of DSD parameters. The path-integrated attenuation (PIA) is also plotted in the figure. The precipitation events occurred from China Standard Time (CST) 00:00 to 08:12 and from CST 10:08 to 12:04, totaling 517 min. The retrieved profile of DSD parameters demonstrates a close correlation with the process of precipitation. According to the profiles of Z, the AGL height of the melting layer increased around the 110th to 120th samples of precipitation, with echo intensity exceeding 35 dBZ, indicating potential convective activity accompanying the precipitation. Based on the statistic results of lighter rain (maximum value of Z under 30 dBZ) and heavier rain (maximum value of Z exceed 30 dBZ) in this day, the retrieved profiles of log ₁₀(N₀) suggest a higher raindrop concentration in the upper levels of precipitation. At top of the precipitation layer, the mean value log ₁₀(N₀) is 5.3 for heavier rain and 4.7 for lighter rain, and the mean of D_m is 0.65 and 0.46 mm. For heavier rain, the vertical variation of the retrieved N₀ and D_m is larger compared to other precipitation samples, the mean absolute value of $(Z_{h_{i+\mathrm{1}}}-Z_h)/\mathrm{\Delta }h$ is 0.01 for log ₁₀(N₀) and 0.002 for D_m. For lighter rain, the value of log ₁₀(N₀) and D_m is 0.009 and 0.001, respectively. Samples 480 to 500 occurred at the end of a brief rainfall in the morning. The maximum value of Z is below 20 dBZ. The retrieved profiles of log ₁₀(N₀) were generally larger than those of other precipitation samples, while the retrieved profiles of D_m indicated smaller raindrop sizes. Based on the calculated PIA and the retrieved D_m profiles, the attenuation is attributed to larger-sized raindrops, with more than 2.6 dB km⁻¹ of attenuation observed in the upper levels of the convective precipitation.

Figure 8 The distribution and the mean values of DSD parameters with height during July to August 2024 at Hongyuan. (a) log ₁₀(N₀). (b) D_m.

4.2 Characteristics of the vertical profiles of DSD parameters

Figure 8 presents the distribution and the mean values of DSD parameters with height during July to August 2024 at Hongyuan. The maximum Z in the profile is used as an indicator of the heavy precipitation. If the indicator is larger than or equal to 30 dBZ, the precipitation sample is identified as heavy precipitation. As shown in Fig. 8a, the mean values of the log ₁₀(N₀) for heavy precipitation are notably higher above 2 km, indicating a larger particle concentration at higher altitudes of heavy precipitation. At altitudes below 1.9 km, the mean values of the log ₁₀(N₀) are lower than those of other lighter precipitation. This is consistent with the prior information that the particle concentration is lower in heavy precipitation. The mean values of log ₁₀(N₀) and D_m for lighter precipitation show vertical variation. According to the statistical results for the D_m parameter, the mean values of D_m for both heavy and lighter precipitation above 2 km are approximately 0.5 mm. This may suggest that raindrops begin to detach from the cloud layer (Barthazy et al., 1998). The mean D_m for heavy precipitation is generally larger than lighter precipitation, and increases as the raindrops fall, reaching a peak at around 0.5 km. This may relate to the equilibrium of coalescence and breakup effects of raindrops during their falling process (Gatlin et al., 2018).

Figure 9 Box plots of the DSD parameters at 210, 500, and 2000 m. (a) N₀. (b) D_m.

Figure 9 presents the diurnal distribution of DSD parameters at heights of 210, 500, and 2000 m. The daytime in the Hongyuan region is longer than the nighttime in summer. The sunrise occurs around 06:20 CST in the morning, and the sunset around 20:30 CST. Usually, the precipitation intensity in this area in the afternoon is higher than that after the cooling at nighttime (Liu and Yin, 2001). The broader distribution for the log ₁₀(N₀) and D_m parameters exhibits around 17:00 CST, likely associated with the convective precipitation that occurs during the afternoon. The mean value of the log ₁₀(N₀) parameter at three different heights between 12:00 and 22:00 CST is lower than other times. This may relate to prior information on the DSD, which suggests that heavy precipitation typically has lower log ₁₀(N₀) values. The D_m have larger mean values between 13:00 and 24:00 CST at three different heights, indicating that the sizes of raindrops are generally smaller during the nighttime. This is consistent with the previously summarized features of DSD of lighter precipitation.

Figure 10 Empirical relations for QPE and attenuation correction using Ka-band radar in the eastern Tibetan Plateau. (a) Z-R relation. (b) D_m-R relation, (c) Z-A relation. (d) D_m-$A/Z$ relation, where A refers to the specific attenuation in dB km⁻¹, Z in mm⁶ m⁻³.

Compared with the lower reaches of the Yangtze River and the precipitation centers in South China, the melting layer of summer rainfall in the eastern Tibetan Plateau is lower, distributed from 0.72 to 2.1 km. The mean value of D_m is approximately 0.6–1.0 mm smaller (Huo et al., 2019; Li et al., 2024). Due to the differences in the features of the DSD parameters in the eastern Tibetan Plateau, empirical relationships for Z-R and D_m-R for QPE using Ka-band radar during summer days and nights are fitted, as shown in Fig. 10a and b. According to the fitted empirical relationships, for the same value of R, the proportion of large raindrops in daytime precipitation is higher. Additionally, the Z-A relationship and the D_m-$A/Z$ relationship that applied single-frequency Ka-band radar to retrieve the precipitation parameters are fitted, as depicted in Fig. 10c and d. The fitted coefficients of the Z-A and the D_m-$A/Z$ relationship in the eastern Tibetan Plateau are compared with the empirical relationships obtained from the OceanRAIN database at all latitudes (Protat et al., 2019). The results indicate that for a measured Z, the specific attenuation is smaller, which is consistent with the smaller raindrop size in the eastern Tibetan Plateau.

5 Conclusions

Using observations from the Ka-band zenith-pointing Doppler radar and the ground-based disdrometer deployed in Hongyuan, the spatiotemporal characteristics of the vertical profiles of DSD parameters in the eastern Tibetan Plateau during summer are studied. A method based on optimal estimation theory is developed for retrieving the vertical profiles of DSD parameters from a Ka-band zenith-pointing radar. The method corrects the attenuation layer by layer, and retrieves the log ₁₀(N₀) and D_m of the precipitation layer from bottom to top. The uncertainty estimation for the retrieval is provided based on the optimal estimation theory. The effectiveness of the optimal estimation algorithms is demonstrated using approximately 5900 synergistic observations from the disdrometer and the Ka-band zenith-pointing radar during the field campaign in the summer of 2024. The average bias of the log ₁₀(N₀) derived from the ground-based disdrometer and the retrieved values from radar observations is 0.12, while the average bias in D_m is −0.01 mm, most of the values observed by the disdrometer falling within the uncertainty range of the retrievals.

According to the retrieved vertical profiles of DSD parameters, the heavy rainfall in this region has higher particle concentration above 2.0 km and a larger particle size at all heights on average. The mean values of D_m are approximately 0.5 mm above 2.0 km. For heavy precipitation, the particle size increases as the raindrops fall, reaching a peak at around 0.5 km. At the bottom of the heavy rainfall, the mean log ₁₀(N₀) is lower than that of lighter precipitation, which is consistent with the prior information derived from the ground-based disdrometers.

In Hongyuan, heavy precipitation primarily occurs between 13:00 and 23:00 CST. During this period, the mean value of log ₁₀(N₀) is lower, and D_m is larger. A broader distribution of the log ₁₀(N₀) and D_m is observed around 17:00 CST, likely associated with the convective precipitation that occurs during the afternoon. The precipitation after nighttime cooling generally has smaller particle sizes and higher particle number concentrations.

To enhance the accuracy of quantitative precipitation estimation and radar retrieval algorithms in the eastern Tibetan Plateau, Z-R and D_m-R relationships suitable for both daytime and nighttime conditions are fitted. In addition, to improving the accuracy of precipitation parameter retrieval using Ka-band radar in the eastern Tibetan Plateau, Z-A and the D_m-$A/Z$ relationship are fitted. The parameter features of the fitted empirical relationships are consistent with the smaller raindrop size observed in the eastern Tibetan Plateau.

Due to the limited information provided by the single-frequency Ka-band radar, this study employed the M-P drop size distribution to avoid the errors introduced by the assumptions and empirical relationships of the GMD in the retrievals. Further evaluation is required to assess the impact of parameterizing the DSD on retrieval accuracy. The DSD characteristic of warm-rain and cold rain is an important topic, and it should be studied in future work. In addition, the retrievals performance and microphysical characteristics under intense convective conditions are not considered, and raindrops are assumed to be spherical in the radar forward model. Further research is necessary regarding retrievals assumptions in these scenarios with Ka-band zenith-pointing radar, such as the treatment of raindrop shape and DSD parameterizations in the retrieval.