Retrieval of Vertical Structure of Raindrop Size Distribution from Equatorial Atmosphere Radar and Boundary Layer Radar

This work develops an algorithm to retrieve the vertical structure of the raindrop size distribution (DSD) of rain from simultaneous observations of 47 MHz Equatorial Atmosphere Radar (EAR) and 1.3 GHz Boundary Layer Radar (BLR) at Koto Tabang, West Sumatra, Indonesia (0.20°S, 100.32°E, 865 m above sea level). EAR is sensitive to the detection of turbulence, and BLR is susceptible to identifying precipitation echo. The EAR Doppler spectrum broadening effects due to turbulence and finite radar beam width were reduced using the convolution process. The Gaussian function was used to model the turbulence Doppler spectrum. A non-linear least-squares fitting method was applied to retrieve DSD parameters. Subsequently, the equations to estimate DSD using this dualfrequency algorithm assume the gamma DSD model to retrieve the distribution from the Doppler spectrum of precipitation echo. The precipitation events on April 23, 2004 on the Coupling Processes in the Equatorial Atmosphere (CPEA-I) project have been analyzed. Results show that the precipitation spectrum obtained using the dual-frequency method is higher, more precise, and wellfitted than the single-frequency method, meaning the dual-frequency method has great potential to be used in observing the microphysical process and remote sensing application analysis of DSD in Indonesia, particularly at Koto Tabang. The analyses show various microphysical processes that occur in the rain, such as coalescence, evaporation, break-up, and condensation. Furthermore, for the purpose of easier remote sensing application analysis of profile DSD characteristics, we use a DSD ΔΖMP parameter. ΔΖMP is a rain rate insensitive DSD parameter representing mean drop size. The trend of ΔZMP is not totally uniform with regards to rain rate and reflectivity factors, with ΔZMP higher in the first half of the event and becoming lower toward the end. This suggests that we have to use different Z-R relations within the event.

of precipitation systems [8]. Due to its wide range of applications, some studies have attempted to measure DSD directly using Joss-disdrometer, two-dimensional video disdrometer/2DVD [9], videosonde and GPS radiosonde [10], precipitation occurrence sensor system [11], and Global Precipitation Measurement (GPM) [12]. Subsequently, the retrieval of the vertical structure of DSD from various atmosphere radar data has been studied by some investigators [13]. They used the simultaneous observation of Very High Frequency (VHF) and L-Band Doppler radar at Shigaraki, Japan. Schafer et al. [14] estimated the DSD from dual-frequency wind profiler spectra using deconvolution and a nonlinear least square fitting technique. In addition, Iwasaki et al. [15] retrieved raindrop and particle size distributions using 14 GHz and 95 GHz radars.
Although there have been several studies on the vertical profile of DSD, research for Indonesia's tropical regions is still very limited. On the other hand, Indonesia receives high rainfall throughout the year. Besides, rainfall and also DSD vary considerably across the region [16]. Unfortunately, the research on the DSD in Indonesia is mostly related to the DSD at ground level [17][18][19][20]. In 2001, the Equatorial Atmosphere Radar (EAR), a VHF band wind profiler, was completed in Kototabang, West Sumatra, Indonesia ( Figure 1). This radar provides an opportunity to observe the vertical structure of DSD in Indonesia. EAR was designed to have the ability to measure the troposphere and lower stratosphere wind field up to 20 km altitude [21]. It is operating at 47 MHz, and capable for detecting two simultaneous signals, one from clear air turbulence and another from hydrometeor particles. Several studies used EAR data or called singlefrequency algorithm to estimate of DSD [18][19][20]. Nevertheless, VHF is limited in retrieval drop sizes smaller than 1 mm in diameter, because the clear air echo has a finite spectral width and dominates the precipitation echo at fall speeds that correlate with small drops [22].

Figure 1. Location of study area
To overcome such a limitation, this work improves the algorithm to resolve the vertical profile of DSD by combining EAR with boundary layer radar (BLR), also known as the dual-frequency algorithm. BLR transmits radio waves with a frequency of 1357.5 MHz, which is a UHF band that is used for observation of wind and hydrometeors in the lower atmosphere, mainly in the planetary boundary layer (PBL). The BLR is more sensitive to receive precipitation echo than that from ambient air motion [23]. Therefore, we use precipitation echo from BLR to estimate DSD. This is expected to provide a more accurate DSD than only using EAR data, as conducted in previous studies [18]. There is a beam broadening difference between EAR and BLR. A correction of beam broadening was carried out to overcome this limitation. Although there have been many studies on dual-frequency algorithms to calculate DSD, such algorithms have never been applied to EAR. In addition, several studies related to dual-frequency algorithms for observing DSD [14,15] used different frequencies from the current study. Thus, this research can be an additional reference for estimating DSD vertical structure using radar.

2-1-1-Equatorial Atmosphere Radar
Equatorial Atmosphere Radar (EAR) is a large-aperture Doppler radar and the first radar in equatorial regions for measuring the troposphere and lower stratosphere wind field up to 20 km altitude and making observations with a resolution of 75-150 m. EAR has a higher sensitivity in providing the information of updrafts and downdrafts of the clear air motions. A detailed description of the EAR is given by Fukao et al. [21].
The data during the Coupling Processes in the Equatorial Atmosphere campaign I (CPEA-I) (10 April-9 May 2004) with the altitude range between 2139 m and 3940 m ASL were analyzed. Due to time-consuming to run the computer program for all rain events, some events were selected to evaluate the performance of dual-frequency algorithm. The selected date is 23 April 2004, when there are simultaneous observations between EAR and BLR (Table 1).
Boundary-Layer Radar (BLR) transmits radio waves with a frequency of 1357.5 MHz, which is used for observation of wind and hydrometeor in the lower atmosphere, mainly in the planetary boundary layer (PBL). The BLR is more sensitivity to receive precipitation echo than the echoes from ambient air motion. Therefore, we use precipitation echo from BLR to estimate DSD. BLR is about 300 m from EAR site. Both instruments are located in Koto Tabang Global Atmosphere Watch (GAW) station, Bukittingi, West Sumatra, Indonesia (0°12′07″ South Latitude -100°19′05″ East Longitude), see Figure 1. More detailed descriptions can be found in several papers [23][24][25].

2-1-2-The Joss-Waldvogel Disdrometer (JD)
Joss-Disdrometer is an instrument which measures DSD on the ground. Joss-Disdrometer was developed by J. Joss and A. Waldvogel and is manufactured by Distrometer Ltd. in Switzerland. It is an impact-type electromechanical counter and able to record the DSD both in very weak and very strong showers (more than 200mm/h) [26]. With diameter ranging from 0.3 to 5.3 mm sorted into 20 classes, it is possible to record up to 200 drops per second [27].

2-1-3-Calibration Data
The radar reflectivity factor (Z), which is measured by BLR, is needed to be recalibrated, by assuming the radar reflectivity factors observed from the Joss-Disdrometer as the reference value. The calibration method used in this study is subtraction between integral of radar reflectivity factor from Joss-Disdrometer and BLR from beginning time (t = a) to ending time (t = b) on one day, calibration value (Cz) as shown by: where Xj is the reflectivity factor of radar (BLR) to be calibrated (mm 6 /m 3 ), Yj is the reflectivity factor estimated from Joss-Disdrometer (mm 6 /m 3 ). Although we use two radars in this dual-frequency algorithm, the most important is to calibrate the BLR data, because reflectivity from BLR will be applied as the initial value in DSD estimation and essential to obtain accurate DSD estimates. The reflectivity factors of BLR from 675 m to 1125 m are averaged for the calibration value. Calibration is performed each day. Calibration value for 23 April 2004 is -0.707dB.

2-2-1-Doppler Spectrum Model
Doppler velocity spectrum of precipitation echo without atmospheric turbulence and wind Sp(v) is expressed by the equation [28]: where C is a constant, N(D) is DSD. The Gaussian function is used to fit the Doppler spectrum of turbulence echo St(v): where p0 is the peak of spectral power, w is the mean vertical wind velocity, σE is the standard deviation EAR spectrum.
In this dual-frequency study, EAR Doppler spectrum of atmospheric turbulence SEAR(v) is given by: where Pn is the noise level on the spectrum, W(v) is the Fourier transform of the auto-correction function of the rectangular time window, and * represents the convolution integral operation. To interface from EAR to BLR estimation in the dual-frequency algorithm, we made a beam broadening correction, which is expressed by: where the standard deviation of the BLR spectrum (σB), horizontal wind speed (vh), EAR beamwidth (θE), BLR beamwidth (θB). The value of vh is retrieved from zonal wind (u) and meridional wind (v) speed.
BLR Doppler spectrum of the raindrop SBLR(v) is expressed as follows:

2-2-2-Non-linear Least Square Fitting and Dual-Frequency Algorithm
Levenberg-Marquardt non-linear least square fitting was used to retrieve DSD parameters. Eight parameters to be estimated by the fitting are four in the EAR estimation turbulence spectrum, i.e., w, σE, p0, and Pn, and four in BLR estimation precipitation spectrum, i.e., m6 (see Equation 10), the slope parameter (Λxy), shape parameter (μ), and Pn.
As the result, we had w and σB from EAR spectrum that useful to interface to the dual-frequency algorithm. Furthermore, on the BLR precipitation spectrum algorithm, μ was obtained by iteration procedure then choose the best μ value with the minimum RMS error between measured and fitted spectra. Finally, the three DSD parameters m6, Λxy, and μ, were estimated.

2-2-3-DSD Model for Dual-Frequency Algorithm
We used the gamma DSD model with three parameters; N0, Λ, and μ as follows: A choice of such parameterization is a similar concept to the "normalized N0" (N0*), proportional to M3/Dm 4 = M3 5 /M4 4 , where Mx is the x-th moment of DSD and Dm is the mass-weighted mean diameter [26]. Mx is given by: where Γ(μ+x+1) is the complete gamma function. Choosing two arbitrary moments, Mx and My, in this study, we used the gamma DSD model, which is expressed by the Equation 10: where my = My/Γ(μ+y+1) and Λxy = (mx/my) 1/(y-x) . In this study, because Λ is dependent on rainrate (R) and radar reflectifity (Z), we use x = 3.67 (moment of R) and y = 6 (moment of Z). DSD parameters are described by my, Λxy, and μ. Λxy indicates the scale parameter Λ obtained from Mx and My. Retrieval of Eq. 4 based on the Doppler spectrum, which is is the terminal velocity. The fitting of the Doppler spectrum can be made effectively by using the DSD parameter, having a high sensitivity to the Doppler spectrum, i.e., m6 [27].

2-2-4-ΔΖMP Parameter
In addition to DSD parameters, we also used a ΔΖMP (dB), as written in Equation 11, to visualize the vertical structure of DSD quantitatively. It is closely correlated with the Z-R relationship, which is applicable in radar remote sensing of rain rate: where R and Z are rain rate and radar reflectivity factor respectively. Since R is approximately proportional to the 3.67 th moment of DSD, M3.67, i.e. R = cR M3.67 [28], so ΔΖMP is given as follows [18]:  (12) where C = -10log10 α, α = 200cR, dBR = 10log10 R and M6 = Z is the 6 th moment of DSD. The meaning of 2. 33 ̅̅̅̅̅̅̅ is the ratio of M6 to M3.67 or the D 3.67 -weighted mean of D 2.33 . ΔΖMP is expected to have similar characteristics as the median drop diameter D0 or mass-weighted mean diameter Dm = M4/M3, except that the positive correlation between D0 or Dm and R is corrected with the term -0.6dBR. Therefore, ΔΖMP has the characteristics; for practical radar application and for a measurement "mean diameter" weighted to intermediate to large drop diameter range since it is derived from M6 and M3.67. ΔΖMP is a rain rate insensitive DSD parameter representing mean drop size [29]. Figure 2 shows the systematic step to estimate DSD using dual frequency radar. Some conclusions for the explanation are also given in sub-chapter 2-2 or study method.

3-1-Comparison Spectrum from Single and Dual Frequency Algorithm
The dual-frequency method shows a well-fitted and adequately higher spectrum for precipitation. The capability of the dual-frequency spectrum to capture rain echo is better than the single frequency as shown in Figure 3. The figure presents the four spectrum comparisons, namely measured dual-frequency, fitted dual-frequency, measured singlefrequency, and fitted single-frequency at heights 2139, 2440, and 2890 m above sea level at 14:08 Local Time (LT) on April 23, 2004. The precipitation spectrum appears well-fitted in the dual-frequency spectra. The fitting results at single and dual frequency obtained a mean correlation coefficient (r) about 0.978 and 0.995, respectively. Results show that the precipitation spectrum obtained using the dual-frequency method is higher, more precise, and well-fitted than the single-frequency method, meaning the dual-frequency method has great potential to be used in observing the microphysical process and remote sensing application analysis of DSD in Indonesia, particularly at Koto Tabang. This is indicated by the Doppler spectrum in Figure 3, which looks higher than the single frequency. At this point, we find that UHF can capture drops smaller than 1 mm, whereas VHF is unable. This is the reason why dual frequency (combined UHF and VHF) will produce a higher precipitation spectrum than single frequency (VHF). The results are in line with those of previous studies [14,22,[30][31][32], where the median absolute error (MAE) between D0 determined by deconvolution methods and D0 computed by the convolution fitting approach is less than 0.15 mm when these dualfrequency retrieval approaches are applied to real measurements recorded during rain episodes [14]. Circles are the single measured-spectrum, blue lines are single fitted-spectrum, plus signs are dual measured-spectrum, and red lines are dual fitted-spectrum.

3-2-Comparison for Several Heights and Times
Analysis of the suitability or good agreement of the results of dual-frequency algorithms is basic to improving microphysical studies. Section 3.1 has shown the better results of the dual-frequency approach. In this sub-section, we compared the DSD from Joss-Disdrometer observations with that obtained from dual-frequency and Marshall-Palmer Conversely, at 13:40 LT, the DSD of the MP model appears narrower than for the Joss-Disdrometer, and occurs until the middle of the event (14:10 LT). This means that different processes caused different DSD characteristics at the ground level, for example, coalescence, break-up, and evaporation. This is in line with the results of the research in Kobayashi et al. [33] where the value of radar reflectivity will increase and decrease in tropical oceans. Radar reflectivity values are influenced by microphysical precipitation factors.    (Table 2). If we assume the mean terminal velocity from the upper level to the ground is 5 m/s, delay time occurs at about 200 s/km. It seems that at the beginning to the middle of the event (13:40-14:00 LT), rain rate increases from 3940 m to 3040 m, then decreases until the final Joss-disdrometer measurement. There should be several processes here. At the beginning of the event (13:40), the DSD becomes slightly broader, and the rain rate increases (3940 m -3040 m), suggesting that condensation occurs [8]. From 3040 m to ground level, DSD is almost constant and rain rate decreases, suggesting that evaporation occurs. The evaporation process is dominant at the beginning of rain events because of surface heat at ground level or the lowertropospheric relative humidity of the equatorial region that makes droplets re-evaporate [33].
From 13:50 to 14:00 LT, since the DSD is broader on the ground than in the upper level, coalescence toward the ground is suggested. From the upper level to middle height (around 3040 m), coalescence and condensation might occur, which makes the rain rate increase. From the middle height to the ground, the evaporation process is dominant, so the rain rate decreases. Marzuki et al. [20] found that the evaporation and updraft correlated with intensive convection while removing small drops from the spectrum in Koto Tabang.
Opposite DSD characteristic profiles are found from the middle to the end of the rain event. More variation processes generate different fluctuations in rain rates from 3940 m to 2290 m; however, the rate increases at ground level. Coalescence and condensation processes may be the main microphysical processes here. In the last event, DSD and rain rate indicate constant values with height, which suggest a saturated precipitation process, i.e., less intense and/or fewer microphysical processes. Such microphysical processes are relevant to the schematic descriptions of the effects of various processes on the shape of the DSD which is described by Rosenfeld and Ulbrich [8].

3-3-Vertical Structure of DSD Parameters and ΔΖMP
Such microphysical processes the dual-frequency algorithm is able to analyze, with DSD profile characteristics over some time ranges and the vertical structure of the DSD from some heights continuously. Figure 5-a shows several DSD parameters, Doppler velocity, spectral width of precipitation echo, reflectivity, rain rate, median drop diameter (D0), and shape parameter (μ), to discuss the microphysical processes which were observed in Koto Tabang on April 23, 2004 from 13:00 -15:00 local time. Figure 5-b plots DSD parameters such as Doppler velocity, spectral width, reflectivity, rain rate, and ΔZMP for the purpose of the remote sensing application analysis.
Microphysical interpretations from measurement parameters give further information about microphysical processes affecting DSD during their fall. Rainfall Doppler velocity is increasing, see Figure 5-a1. Similarly, spectral width and radar reflectivity are increasing. We can infer that a melting process occurs here, where ice crystals turn to droplets of rain. The highest radar reflectivity factor occurs in the middle of the rain event. In Figure 5-a1, 5-a2, and 5-a3 describe microphysical rain processes from 1165 m to 9750 m ASL from BLR data. Dual-frequency algorithm output parameters are rain rate in dBR, D0 and μ, see Figures 5-a4, 5-a5 and 5-a6, respectively, with a height range from 2000 m -4000 m ASL. For comparison, we also show the Joss-disdrometer rain rate in Figure 5-a7 at ground level.
Figures 5-a1 to 5-a7 are layers to describe the microphysical rain processes from upper to ground level. The atmosphere layer has a freezing level occurring at 0 o C height. The freezing level is dominated by ice crystals which play an important role in the precipitation. Going down about 4 km, there is a microphysical process of melting. From 4 km, we can analyze microphysical processes of the DSD because the precipitation changes to rain droplets. To analyze how microphysical rain processes occur from melting level to the ground, we need to understand the terminal velocity which correlates with the delay time of droplets to reach the ground. As a brief explanation, one of the most common v0-Z relationships for rain is obtained by v0 = 3.84 Z 0.0714 , where v0 is the mean terminal velocity of rain (m/s) and Z is the radar reflectivity factor (mm 6 /m 3 ) [34]. For this data (April 23, 2004), after calculating Z values that include heavy rain, we assume v0 is about 5 m/s, thus the delay time is about 6 minutes. For ease of understanding the analyses, we have placed lines in Figure 5 to explain the microphysical processes that occur in the event. Note that when using the word "line", we refer to a time range of droplets falling down.
In figures 5-a4 and 5-a5, line 1 (13:30-13:38) and line 2 (13:38-13:46) show that rain rate is almost constant along these lines, and D0 seems constant in line 1 and increases in line 2. The DSDs along line 2 are broader than line 1, and describe coalescence and condensation. Condensation still continues, so the rain rates increase at the beginning of line 3. In the middle of line 3, rain rates and D0 decrease, showing that evaporation processes occurs. In the line 4, D0 and rain rates appear constant, but D0 is broader than in line 3, meaning that coalescence occurred from line 3 to 4. D0 and rain rates decrease again in line 5, suggesting that break-up and evaporation occurs here. At the end of the event, lines 6, 7, and 8 show D0 and rain rates of constant value, indicating less intense microphysical processes because of the saturated processes of rain. Based on the fact that IRPs for radar remote sensing and microwave communication applications are mainly proportional to the 3 rd -6 th moments of DSD [7], Figures 5-b can be used for remote sensing application. Considering that the most direct property relating to DSD for radar remote sensing is the Z-R relation, we examine Z-R relation with the parameter ΔZMP. Figures 5(b1)-(b3) show Doppler velocity, spectral width, and reflectivity factor from BLR data. Z and R for the rain region (2-4 km ASL) are shown in Figures 5-b4 and 5-b5, and it appears that the contours of Z and R are similar. The highest value is in the middle (14:30-15:00), second is at the beginning moment (13:20-13:40) and the lowest at the end of the event. Figure 5-b6 shows a ΔZMP plot which seems to have small values of DSD about 0 to -15 dB. Noting that positive and negative values of ΔZMP indicate broad and narrow DSDs, respectively, in comparison with a DSD generating Z-R relationship (Z = 200R 1.6 ). The trend of ΔZMP is not totally the same as with rain rate and reflectivity factor. ΔZMP is higher in the first half of the event and becomes lower toward in the end. This suggests that we have to use different Z-R relations within the event.

4-Conclusion
The results of this study show that the dual-frequency method of combining turbulence echo from the EAR with precipitation echo from the BLR provides an accurate DSD estimation. The fitting results at single and dual frequency obtained a mean correlation coefficient (r) of about 0.978 and 0.995, respectively. Results show that the precipitation spectrum obtained using the dual-frequency method is higher, more precise, and well-fitted than the single-frequency method. This can also be seen from the accuracy of the turbulence and precipitation spectrum fittings and the comparison of DSD values with those of the Joss-Waldvogel Disdrometer. It means the dual-frequency method has great potential to be used in observing the microphysical process and remote sensing application analysis of DSD in Indonesia, particularly at Koto Tabang. This can be seen from the evaluation of rain on April 23rd, 2004. Several retrieved DSD parameters give further information about possible microphysical processes affecting DSD during its fall. The evolution of the DSD in the rain column can be clearly observed. Various microphysical processes in the rain, such as coalescence, evaporation, break-up, and condensation, can be seen from the retrieved DSD vertical structure. This is in line with the characteristics of ΔZMP. Thus, the dual-frequency algorithm can be used to observe the vertical structure of DSD based on EAR and BLR. This work is limited to a rain event. The implementation of the dual-frequency method to characterize the vertical structure of DSD in Kototabang for long data records of the EAR is being carried out, and the results will be published in a subsequent paper.

5-2-Data Availability Statement
The data presented in this study are available on request from the corresponding author.

5-3-Funding
This study was supported by 2021 World Class Research Grants from Ministry of Research, Technology and Higher Education/Ministry of Education, Culture, Research, and Technology (Contract no: 104/SP2H/LT/DRPM/2021). The observation during the CPEA project has been supported by Grant-Aid for Scientific Research on Priority Areas funded by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT).

5-4-Acknowledgements
Thanks to Research Institute for Sustainable Humanosphere (RISH), Kyoto University and the National institute of Aeronautics and Space of Indonesia (LAPAN) for EAR and BLR data. Special thanks to the late Prof. Toshiaki Kozu for fruitful discussion.

5-5-Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this manuscript. In addition, the ethical issues, including plagiarism, informed consent, misconduct, data fabrication and/or falsification, double publication and/or submission, and redundancies have been completely observed by the authors.