Effect of Vertical Air Motion on Disdrometer Derived Z - R Coefﬁcients

: For synoptic-scale motions the vertical velocity component is typically of the order of a few centimeters per second. In general, the vertical velocity is not measured directly but must be inferred from other meteorological ﬁelds that are measured directly. In the present study, a Joss–Waldvogel disdrometer was used in order to establish the drop size distributions (DSD) at Athalassa, Cyprus. Data from a radiosonde station co-located with the disdrometer were also collected which were subsequently used to derive estimates of vertical velocities. Meteorological ﬁelds, including vertical velocities, were extracted from an atmospheric reanalysis, for an area centered over the disdrometer and radiosonde station instrumentation. The disdrometer data were used to determine the Z-R disdrometer derived coefﬁcients, A and b , where Z = A R b . To model the vertical air effect on the Z-R disdrometer derived coefﬁcients an idealistic notion of ﬂux conservation of the DSD is adopted. This adjusted DSD (FCM-DSD) is based on the exponential DSD and is modiﬁed by the relationship between drop terminal velocity ( D ) and vertical air speed w . The FCM-DSD has a similar appearance to the popular gamma DSD for w < 0. A clear segregation is seen in the A - w plane for both data and model. The data points are also clearly segregated in the b - w plane, but the model points are on opposite sides of the w = 0 line. It is also demonstrated that vertical velocities can be extracted from radiosonde data if initial balloon volume is accurately measured, along with an accurate measurement of the mass of the complete radiosonde-balloon system. To accomplish this, vertical velocities from radiosonde data were compared to reanalysis vertical velocity ﬁelds. The resulting values of initial balloon volume are found to be within the range of measured values.


Introduction
One meteorological variable which has a profound effect on the weather is the vertical velocity of the atmospheric air (hereafter denoted by w, in m s −1 ) [1][2][3]. Through sustained vertical motions, rising moist air cools adiabatically forming precipitating cloud [4], sometimes accompanied with more notable phenomena, like lightning and thunder [1]. Rising motions can lead to a steeper environmental temperature lapse rate [5] and strengthening of cyclonic systems [1,3]. On the contrary, descending air is heated by diabatic compression leading to cloud dissipation, damping of precipitation, clearer skies, and finer weather.
Bearing in mind the importance of vertical motions in the atmosphere, this study presents the methodology and results of an attempt to deduce information about the character of the Despite their deficiencies, both methods are widely applied adopting an isobaric coordinate system, (i.e., a coordinate system in which the vertical coordinate is atmospheric pressure instead of geometric height) so that ω = dp/dt is inferred rather than w = dz/dt (p denotes atmospheric pressure, z the geometric height above the earth's surface, and t is time). In an isobaric coordinate system where as explained above pressure replaces height as the vertical coordinate, the pressure tendency, namely ω, plays the role of vertical velocity w = dz/dt in the Cartesian coordinate system; hence, ω = dp/dt is often referred to as the ω-vertical velocity following the motion.
Most commonly, vertical velocities available in data repositories of atmospheric variables are given in terms of ω. For example, vertical velocities in the ERA-Interim database (employed hereafter) maintained by the European Centre for Medium Range Forecasts (ECMWF) are given as pressure tendencies, ω = dp/dt. Nevertheless, as shown later in this paper, converting from ω-velocity to w-velocity is feasible.
An alternative approach for estimating ω fields in the atmosphere, based on the so-called omega equation [3], that is free from the difficulties inherent in the above mentioned kinematic and adiabatic methods may be adopted. Nevertheless, this approach is also not free from deficiencies, despite its wide utilization [9].
The problem of direct observability of vertical velocities in the atmosphere is quite difficult and despite efforts to establish suitable directly observing systems, the problem is far from reaching a globally satisfactory level of coverage. On the one hand, lidar, sodar and wind profiling radar are capable of measuring vertical velocities (e.g., [10,11]), but only under conditions of presence of cloud droplets or ice crystals or aerosol particles, whereas under clear air conditions, no vertical velocities can be determined [12]. Hence, both lidars and cloud radars are practically used to derive information about particle velocities. On the other hand, wind profilers can measure vertical velocities in clear air [13]. However, although they have been explored for this purpose [14], their global-scale distribution is still very limited. Instruments like the above that can be used to observe directly vertical velocities of the air, are still under exploration mostly in campaign experiments [12].
SODAR (SOnic Detection And Ranging) comprise another technique that is used to measure wind speed at various heights above the ground; this wind profiler is based on the scattering of sound waves by atmospheric turbulence. SODAR has also been around for several decades [15].
The application of proxy techniques in estimating vertical velocities is also under investigation. For example, trace gas observations from satellite remote sensing instruments are compared to modeled trace gas distributions [16].
In any case, no global observing network for vertical velocities of the atmospheric air is in place and it is not foreseen to be established any time soon. This gap leaves ample room for further research focusing on the vertical velocity issue. The continual scientific interest in the determination of the vertical velocity fields is reflected in recent studies, such the one by Stepanyuk [17] where the authors present a comparison between ω derived from the ERA-Interim database (discussed immediately below) and values of ω derived from the generalized omega relationship.
Reanalyses of meteorological fields comprising long-term data sets with high spatiotemporal resolution are suitable for climate studies (e.g., [18]); the production of homogeneous data sets is at the core of such endeavors. However, comparing five reanalyses, Iwasaki et al. [19] noted important differences, even in zonally averaged vertical velocities. Nevertheless, it is worth noting that several aspects in the latest generation reanalyses have been enhanced [20]. Despite the improvements, vertical velocities in reanalyses still suffer from inaccuracies, especially in relation to particular applications [18].
Employing disdrometer DSD data to estimate the vertical air velocity w has been approached in different ways, depending on the nature of the disdrometer. The method discussed in this current work is based on a procedure presented in previous work [21], using impact disdrometer data. In this case, the disdrometer extracted A-b coefficients, where Z = AR b , seem to be influenced by the magnitude of w and especially the sign. In a more recent work, Kim and Song [22] have compared w measured by a three-axis ultrasonic anemometer and by a collocated laser disdrometer. The w considered by Kim and Song [22] is not the same w considered by Lane et al. [21] and in this current paper. In the present work, w is an average value from the surface to a reference height of z ref = 10 km. In the work by Kim and Song [22], w is the value measured at the surface. On a flat surface, w should be zero, but because Kim and Song studied areas near mountain slopes, the vertical component of air motion at the surface can be non-zero.
The aim of the present research comprises an attempt to deduce information about the character of the vertical air motion w from the disdrometer data. In this respect, the results of combining measurements from a Joss-Waldvogel disdrometer, radiosonde measurements and fields from atmospheric reanalysis are presented and discussed. The present paper aspires to contribute towards the better understanding of the long-standing vertical velocity estimation issue through the combined interplay of these three sources of data.
Section 2 presents the data used in this analysis together with the preprocessing approaches adopted. The methodology adopted is outlined in Section 3, whereas, the research outcomes are presented in Section 4 and discussed in Section 5. The paper concludes with Section 6. A list with the abbreviations and symbols used in this paper is given in Appendix A.

Data and Their Preprocessing
In this study, three different sources were used: disdrometer data, co-located radiosonde data, and data from the European Centre for Medium Range Forecasts (ECMWF). These data sources are discussed in the following together with the preprocessing procedures adopted in order to formulate them for further use in this research.
This study sets the theoretical foundations in reaching the aim of this research and only a limited number of cases were investigated. Eleven days in the period 2011-2014 were selected to showcase the aim of this research. The dates selected are the same used in a previous study by Lane et al. [21] and are shown in Table 1 of the present paper. In searching for candidate dates, the authors have tried to balance the updraft and downdraft cases and a variety of DSD types was also pursued. The relation to stratiform-convective precipitation was described in the paper by Lane et al. [21] in which all vertical speeds were simulated values and parameter fits. In the current work, matching ERA-Interim reanalysis data was used to get vertical velocity data, together with collocated radiosonde data.

Disdrometer Data
The disdrometer data were recorded by using a Joss-Waldvogel impact disdrometer located on the roof of a building at Athalassa, Cyprus (35. convert impulse size to drop diameter. The Joss-Waldvogel impact disdrometer is able to record drop diameters from 0.3 mm to 5.5 mm in 10-second intervals, thereby allowing for the establishment of the drop size distribution (DSD) representing this range of drop sizes [23,24]. As long as the manufacturer's calibration certificate has not expired, the data is considered to be reliable.

Radiosonde Data
The data from the radiosonde station at Athalassa (next to the disdrometer) are collected for each of the days in the study. The equipment used is a Vaisala DigiCora ® Sounding System MW41 which is GPS-based. The data are automatically collected by the ground station and are registered for every two seconds, as the radiosonde ascends freely in the atmosphere. A sample of a file containing radiosonde data is given in Appendix B. These data consist of the pressure (p), geopotential height (Z h , in geopotential meters, gpm), temperature (T, in • C), relative humidity (RH, in %), Dew-point temperature depression (i.e., the difference between temperature and dew-point temperature, in • C), virtual temperature (in • C), temperature lapse rate (in • C km −1 ), and balloon ascent rate (w r in m s −1 ).
Vertical air motion is deduced from the radiosonde data by comparing its GPS-based ascent velocity to the theoretical ascent velocity of a balloon in still air. Wang et al. [25] developed a procedure which results in Equation (1), yielding a calculated balloon still air ascent velocity, w c V 0 is the initial balloon volume at the surface before release, ρ 0 is the air density at the surface, ρ is air density as a function of height, m T is the mass of the complete radiosonde-balloon system (m T = 0.9 kg, for the radiosonde system considered in this work), and g 0 is the global average of the acceleration due to gravity at MSL. The radiosonde reports several variables during its ascent as a function of time, which can easily be transformed into a function of height. Air density is neither measured by the radiosonde sensors nor reported in the radiosonde data file; however, density can be calculated from the variables that are being reported in the radiosonde data file, as explained immediately below. From the ideal gas law, air density can be expressed as a function of pressure (p), temperature (T in K) and the gas constant for moist air (R m ) where, R m is defined by with R d (287.058 J kg −1 K −1 ) and R v (461.5 J kg −1 K −1 ), denoting the gas constants for dry air and water vapor, respectively; q denotes the specific humidity (in kg kg −1 ). In the radiosonde data, however, specific humidity is not reported; instead the relative humidity (i.e., RH) is reported. To convert from relative humidity to specific humidity, we use the following relationship which is based on the Clausius-Clapeyron equation where, T 0 is a reference temperature (typically taken to be 273.16 K). Hence, specific humidity is Atmosphere 2019, 10, 77

of 21
These values provide the preprocessing needed to compute Reynold's number Re, coefficient of drag C D , and calculated still air vertical ascent speed w c . These three values are coupled and require a numerical algorithm to compute, as discussed below.
Ultimately, Equation (1) is defined for our needs as a function of height z and it should work for any planet with an atmosphere. If there is no atmosphere, Equation (1) fails to predict anything useful, since it is really the solution to a buoyancy problem such that V 0 ρ 0 > m T . Equation (1) seems straightforward, but there is one problem: the coefficient of drag C D is a function of Reynold's number, which is a function of w c . Wang et al. [25] take an easy way out and define C D to be a constant during the entire ascent. The value they use is in the range of 0.4 to 0.5. Sóbester et al. [26] go one step further and define C D to be a simplified function of Reynold's number where, C 1 = 0.425, C 2 = 0.225, Re 1 = 329600, and Re 2 = 365900. This improvement comes at a cost, since Re is a function of w c and C D where, ρ and w c are from Equation (1), D b is the balloon diameter (which can be calculated directly from balloon volume V), and µ is the dynamic viscosity of the air (in kg m −1 ·s −1 ) given by the formula where, the constants are: A n = 1.862 × 10 −7 and β = 0.8062. Equations (1), (6), and (7) represent a coupled set of equations that must be solved iteratively. A simple recursive solution to the above set of equations is to define an error ε n which is then used to recursively find a solution to Re, C D , and w c . The iterative algorithm to find a solution starts by calculating C Dn from Re n−1 , using Equation (6). Next, w c n is calculated using C Dn and Equation (1). An error is calculated which is then used to update Re Re n = Re n−1 − γε n (10) where, γ is the recursion gain factor set equal to 0.4 in this work. Equations (9) and (10) are computed N times with an initial guess for Re 0 = 10 5 . After N = 100 iterations, the value of ε n drops to < 10 −4 percent for all cases tried, resulting in the convergence of vertical speed at radiosonde sample time.

ECMWF ERA-Interim Reanalysis Data
In this paper, the ERA-Interim reanalysis dataset was utilized. The data is provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) via its internet portal and it will be shown below how the meteorological fields retrieved are used in this work to estimate the vertical velocity, denoted by w (in m s −1 ).
ERA-Interim is a global atmospheric reanalysis from 1979, continuously updated in real time. The ERA-Interim atmospheric model and reanalysis system uses cycle 31r2 of ECMWF's Integrated Forecast System (IFS) [27]. The ERA-Interim reanalysis is described by Dee et al. (2011).
The ERA-Interim analyses produced and distributed by ECMWF, refer to an atmospheric model based on a hybrid vertical coordinate system [27,28]. Such a hybrid system consists of a terrain following coordinate σ = p/p s near the surface (where, p is pressure and p s is the surface pressure), Atmosphere 2019, 10, 77 7 of 21 but with a gradual transition to a pressure coordinate with height. In addition, archived ERA-Interim products comprise reanalysis fields from the atmospheric model evaluated on standard pressure levels; this latter form of data have been retrieved from ECMWF's portal and are subsequently utilized in the present work.
The retrieved data refer to the finest grid available for the ERA-Interim, namely 0.125 • × 0.125 • (which over Cyprus corresponds to a grid with approximate distances along latitude and longitude of 13.9 km and 11.4 km, respectively). The data cover the area from 30 to 40 • N and from 30 to 40 • E. In the vertical, data for the 23 standard pressure levels in the range from 1000 to 200 hPa.
For the needs of the study here, a 1 • × 1 • sub-area centered at the disdrometer site at Athalassa has been considered, as shown in Figure 1. This area is bounded by meridians 33.0 • E and 34.0 • E and by latitude circles 34.6 • N and 35.6 • N. Because of the 0.125 • × 0.125 • ECMWF grid spacing, the sub-area that is extracted from the dataset needs to be no smaller than about 0.2 • × 0.2 • , otherwise the number of data points becomes too sparse for most calculations. With a 0.2 • × 0.2 • region, at least one data grid point will exist for every 1000 m increment above the surface.
Atmosphere 2018, 9, x FOR PEER REVIEW 7 of 21 meteorological variables, including vertical air speed. The following meteorological variables were retrieved from the ERA-Interim database for each grid point and each pressure level: eastward and northward wind components (u, and v, respectively, in m s −1 ), temperature (T, in K), geopotential (Φ,in m −2 s −2 ), specific humidity (q, in kg kg −1 ), and vertical velocity (ω = dp/dt, in Pa s −1 , where p is pressure and t is time; as explained above, ω is defined as pressure tendency in a system where pressure is considered as the vertical coordinate). The vertical velocity is calculated from the ERA-Interim data using the following relationship [3,29] w g where, ρ is density (in kg m −3 ) and g is the acceleration due to gravity. Substituting ρ from Equation (2), the equation used for calculating the vertical velocity becomes where, m R is defined by Equation (3).
The acceleration due to gravity is adjusted for latitude, φ, and height above the MSL, z, as follows: (https://www.sensorsone.com/local-gravity-calculator/) where, IGF is the International Gravity Formula (in m s −2 ), given by the relationship [30,31] and FAC is the Free Air Correction (in m s −2 ), given by the relationship From the definition of the geopotential Φ(z) at height z as the work required to raise a unit mass to height z from MSL, and from the conversion of the geopotential into geopotential height ECMWF provides daily atmospheric data corresponding to a horizontal latitude-longitude grid and are available for every 6 h. This vast data set provides the spatiotemporal fields of many meteorological variables, including vertical air speed. The following meteorological variables were retrieved from the ERA-Interim database for each grid point and each pressure level: eastward and northward wind components (u, and v, respectively, in m s −1 ), temperature (T, in K), geopotential (Φ, in m −2 s −2 ), specific humidity (q, in kg kg −1 ), and vertical velocity (ω = dp/dt, in Pa s −1 , where p is pressure and t is time; as explained above, ω is defined as pressure tendency in a system where pressure is considered as the vertical coordinate).
The vertical velocity is calculated from the ERA-Interim data using the following relationship [3,29] w = − ω ρg (11) where, ρ is density (in kg m −3 ) and g is the acceleration due to gravity. Substituting ρ from Equation (2), the equation used for calculating the vertical velocity becomes where, R m is defined by Equation (3). The acceleration due to gravity is adjusted for latitude, ϕ, and height above the MSL, z, as follows: (https://www.sensorsone.com/local-gravity-calculator/) g = g(ϕ, z) = IGF + FAC (13) where, IGF is the International Gravity Formula (in m s −2 ), given by the relationship [30,31] and FAC is the Free Air Correction (in m s −2 ), given by the relationship From the definition of the geopotential Φ(z) at height z as the work required to raise a unit mass to height z from MSL, and from the conversion of the geopotential into geopotential height where g 0 is the global average of acceleration due to gravity at MSL, g 0 ≡ 9.80665 m s −2 , the geopotential height Z h is derived. In the troposphere and lower stratosphere, the geopotential height Z h is numerically almost identical to the geometric height z, therefore the former can be used in Equation (15) to calculate the free air correction. Figure 2 shows 3D contour plots of 1 • × 1 • ECMWF vertical wind data for each of the radiosonde observation times: 06:00 UTC and12:00 UTC 17 April 2013. The legend colors correspond to the range of values of w. One notable characteristic of these plots is the apparent uniformity of air speeds at various height levels. This uniformity removes any significant dependence on region's window size, as long as a minimum window of 0.2 • × 0.2 • is used. A maximum window of 1 • × 1 • is likely to provide minimum variation at the horizontal levels. All of the variation is in the vertical direction. This appears to be true at the 1 deg × 1 deg scale and smaller. Note that because of the spacing of ECMWF grid points, microscale features are not likely to show up. where 0 g is the global average of acceleration due to gravity at MSL, 0 g ≡ 9.80665 m s −2 , the geopotential height h Z is derived. In the troposphere and lower stratosphere, the geopotential height h Z is numerically almost identical to the geometric height z, therefore the former can be used in Equation (15) to calculate the free air correction. Figure 2 shows 3D contour plots of 1° × 1° ECMWF vertical wind data for each of the radiosonde observation times: 06:00 UTC and12:00 UTC 17 April 2013. The legend colors correspond to the range of values of . One notable characteristic of these plots is the apparent uniformity of air speeds at various height levels. This uniformity removes any significant dependence on region's window size, as long as a minimum window of 0.2° × 0.2° is used. A maximum window of 1° × 1° is likely to provide minimum variation at the horizontal levels. All of the variation is in the vertical direction. This appears to be true at the 1 deg × 1 deg scale and smaller. Note that because of the spacing of ECMWF grid points, microscale features are not likely to show up. In order to extract a useful and meaningful value of from the ECMWF data, all data at 1000 m levels are averaged to create wind profiles, as shown in Figure 3. The wind profile is integrated from the surface to form an accumulated average of the wind data profile. Finally, for comparison to the disdrometer derived A and b parameters as discussed below, the accumulated wind profile at a reference level is used. The reference level chosen for this comparison is zref = 10,000 m. The zref seems to provide a consistent comparison to the disdrometer A-b data; (zref) is simply the average of the wind profile from z = 0 to zref. Table 1 summarizes these results for 11 days of disdrometer-ECMWF data. In order to extract a useful and meaningful value of w from the ECMWF data, all data at 1000 m levels are averaged to create wind profiles, as shown in Figure 3. The wind profile is integrated from the surface to form an accumulated average of the wind data profile. Finally, for comparison to the disdrometer derived A and b parameters as discussed below, the accumulated wind profile at a reference level is used. The reference level chosen for this comparison is z ref = 10,000 m. The z ref seems to provide a consistent comparison to the disdrometer A-b data; w(z ref ) is simply the average of the wind profile from z = 0 to z ref . Table 1 summarizes these results for 11 days of disdrometer-ECMWF data.
In order to extract a useful and meaningful value of from the ECMWF data, all data at 1000 m levels are averaged to create wind profiles, as shown in Figure 3. The wind profile is integrated from the surface to form an accumulated average of the wind data profile. Finally, for comparison to the disdrometer derived A and b parameters as discussed below, the accumulated wind profile at a reference level is used. The reference level chosen for this comparison is zref = 10,000 m. The zref seems to provide a consistent comparison to the disdrometer A-b data; (zref) is simply the average of the wind profile from z = 0 to zref. Table 1 summarizes these results for 11 days of disdrometer-ECMWF data.

DSD Flux Conservation Model
An approach to modeling the vertical air effect on the Z-R disdrometer derived coefficients is to start with the very idealistic notion of flux conservation of the DSD, where the DSD flux, F(D), is defined as  (27) and (28) in Lane et al. [21]). The DSD measured by a disdrometer at the surface is N (D) with flux v (D)N (D), which can be found by equating the flux at the surface and F(D) in the vertical air column

Drop Velocity
Equation (18) is an oversimplification in more than one ways, but it provides a starting point for calculating the effect of vertical air motion on the surface DSD. The most glaring problem with Equation (18) as this is formulated, in addition to the problem of requiring ideal non-interacting drops, is that the right-hand side of this equation turns negative when w is positive (updraft) for some drop sizes below a value of D 0 (w). This lower limit of D can be found as the solution to v(D 0 ) = w. We can use any approximation of v(D), including a power law, to find a D 0 and demonstrate the dependence of the coefficients A and b (based on the Z-R relationship, as explained below) on w.
We will use the following approximation for v (D) where, a 1 = 4.67, a 2 = −0.789, a 3 = 0.0441 (see [21]). The w 0 and λ are model parameters that attempt to approximate the behavior of drops that are approaching the surface with velocities different than terminal velocity due to the finite distance needed to accelerate/decelerate from the vertical air column. The parameter w 0 is in most cases set equal to w. The last term of Equation (19) may be especially useful in the case of a laser disdrometer derived Z-R since v (D) is measured directly. Since the disdrometer at Athalassa, Cyprus, is an impact disdrometer, we will not attempt to utilize the last term, which can be done by setting w 0 = 0 (or setting λ to a very large value). With the definition of drop velocity given by Equation (19), the still air terminal velocity v(D) = v (D) for w 0 = 0. Solving for D 0 , where v(D 0 ) − w = 0 from the numerator of Equation (18), can be accomplished with some effort using standard root solving strategies (or with minor effort using Mathematica ® ). The exact solution was found using Mathematica. It is a long, unwieldy result. To simplify the exact Mathematica solution, a set of data was created from the exact formula, then that data was fitted to a third order polynomial where, α 1 = 0.2502, α 2 = −0.01324, α 3 = 0.003405, and H(w) is the unit step function (also referred to as the Heaviside step function). The step function is required to prevent negative drop diameters and is defined as D 1 and D 2 are defined as the lower and upper limits of the disdrometer sensitivity. Using the Joss-Waldvogel disdrometer limits of D 1 = 0.3 mm and D 2 = 5.5 mm [32], Equation (19) leads to the plot of Figure 5.

DSD Flux Conservation Model
An approach to modeling the vertical air effect on the Z-R disdrometer derived coefficients is to start with the very idealistic notion of flux conservation of the DSD, where the DSD flux, ( ), is defined as where, ( ) is drop terminal velocity and ( ) is the DSD. In this ideal world of flux conservation, drops do not interact, and they maintain their size and population during their lifetime and subsequent fall to the surface. The only characteristic that can change in time or location is drop velocity. Since ( ) is conserved and is therefore a constant for any , any change in ( ) must be accompanied by an inverse change in ( ).  To understand how flux conservation model (FCM) will lead to a change in the DSD spectra, refer to Figure 4, where a DSD described by ( ) is falling with a combined terminal velocity ( ) and vertical air motion (see also Equations (27) and (28) in Lane et al. [21]). The DSD measured by A severe shortcoming of Equation (18) is the implication that no drops will exist for the disdrometer to measure below D 0 (w). However, that is seldom the case (one might say never the case), since disdrometers are not likely to have a completely depleted distribution for all small drops below some drop size. A modification can be made to Equation (18) to account for drop breakup that partially repopulates the size distribution below D 0 (w) in the region above the surface, but below the height where w is constant where, x ≡ D − D 0 (w). The prime has been dropped from the LHS of Equation (22) since, from this point forward, N(D) will correspond to the theoretical DSD at the disdrometer. where, ≡ − ( ). The prime has been dropped from the LHS of Equation (22) since, from this point forward, ( ) will correspond to the theoretical DSD at the disdrometer. Combining Equations (19) and (22) results in what will henceforth be referred to as the FCM-DSD, which is the theoretical DSD experienced by the disdrometer, based on the FCM. Figure 6a shows the FCM-DSD of Equation (22)

Z and R as DSD Moments
Computing Z and R from the DSD model of Equation (22) is nothing more than an exercise in computing integer moments of the DSD (see also [34]) Combining Equations (19) and (22) results in what will henceforth be referred to as the FCM-DSD, which is the theoretical DSD experienced by the disdrometer, based on the FCM. Figure 6a shows the FCM-DSD of Equation (22) with N 0 = 8000 (mm −1 m −3 ), β = 0.9, w 0 = 0, with three values of w and two values of Λ, i.e., the Marshall-Palmer rainfall rate parameter [33]. Figure 6b  where, ≡ − ( ). The prime has been dropped from the LHS of Equation (22) since, from this point forward, ( ) will correspond to the theoretical DSD at the disdrometer. Combining Equations (19) and (22) results in what will henceforth be referred to as the FCM-DSD, which is the theoretical DSD experienced by the disdrometer, based on the FCM. Figure 6a shows the FCM-DSD of Equation (22)

Z and R as DSD Moments
Computing Z and R from the DSD model of Equation (22) is nothing more than an exercise in computing integer moments of the DSD (see also [34])

Z and R as DSD Moments
Computing Z and R from the DSD model of Equation (22) is nothing more than an exercise in computing integer moments of the DSD (see also [34]) Note that the limits of the moment integral are not shown as 0 to ∞, but are defined as the Joss disdrometer detection limits D 1 and D 2 ; Z is then computed directly from the 6th moment: Z = M 6 . R is an integral of the product of drop velocity v(D), drop volume, and N(D). If Equation (19) is used with w 0 = 0, then R becomes the sum of the fourth, fifth, and sixth moments, weighted by the coefficients in Equation (19 where, for k = 1,2,3, b k = a R a k , and a R = 0.0036 π 6 which sets R in standard units of mm where, Γ I is the incomplete gamma function, as a result of the finite limits on the integral in Equation (23). Unfortunately, the form of the DSD based on Equation (25) does not lead to convenient forms for M n , except in the special case of w = 0. In general, it is necessary to numerically evaluate the moments to attain values for Z and R.

A-b Coefficients
A parametric plot of Z versus R over a range of Λ with a fixed N 0 generates nearly straight lines on a log-log plot, (as discussed below in conjunction with Figure 7). This inherent property of nearly straight lines is the source and motivation for expressing Atmosphere 2018, 9, x FOR PEER REVIEW 12 of 21 Note that the limits of the moment integral are not shown as 0 to ∞, but are defined as the Joss disdrometer detection limits D1 and D2; Z is then computed directly from the 6 th moment: Z = . R is an integral of the product of drop velocity ( ), drop volume, and ( ). If Equation (19) is used with = 0, then R becomes the sum of the fourth, fifth, and sixth moments, weighted by the coefficients in Equation (19) = + + (24) where, for k = 1,2,3, = , and = 0.0036 which sets R in standard units of mm h −1 . If ( ) is an exponential distribution where, is the incomplete gamma function, as a result of the finite limits on the integral in Equation (23). Unfortunately, the form of the DSD based on Equation (25) does not lead to convenient forms for , except in the special case of = 0. In general, it is necessary to numerically evaluate the moments to attain values for Z and R.

A-b Coefficients
A parametric plot of Z versus R over a range of with a fixed generates nearly straight lines on a log-log plot, (as discussed below in conjunction with Figure 7). This inherent property of nearly straight lines is the source and motivation for expressing This can be expressed as log = log + log (28) This can be expressed as The coefficient set {logA, b} is found from the linear fit of logR versus logZ. Equation (27) is an approximation that may be good under a limited set of circumstances. This conclusion becomes obvious by inspection of the DSD moment M n in Equation (26). Two requirements would need to be satisfied for Equation (27) to be exact: N 0 would need to be a constant for all time and drop terminal velocity v(D) would need to be accurately represented by a power law. Neither of these requirements are satisfied, except under very limited circumstances. Nevertheless, Equation (27) can be very useful as long as it is understood that it is just an approximation. Figure 7 shows the Z-R curves generated with a parametric plot by running Λ through a range of values (approximately 1 through 7). Three curves are shown from the FCM-DSD model with three values of w, N 0 = 8000 (mm −1 m −3 ), β = 1, and w 0 = 0. The limits of integration D 1 (w) and D 2 as shown in Figure 5, are D 1 (0) = 0.3 (mm) and D 2 = 5.5 (mm). In addition, the National Weather Service (NWS, United States of America) standard Z-R (see [35]) is shown for comparison. Figure 8 shows two sets of curves for comparison with widely varying values of N 0 . From these figures it is obvious that the Z-R line shifts with values of w, corresponding to a different value of A. There are also accompanying changes in the slope. Changing N 0 and β also shift the Z-R lines, producing different values of A and b. Setting β = 0 turns off the effect w has on the FCM-DSD (see Equation (22)). The short dotted and long dotted lines will converge to the solid line for each N 0 shown in Figures 7 and 8. This parameter has the effect, on the DSD model of Equation (22), of mixing in a percentage of the standard exponential DSD from 0 to 100%. The coefficient set {logA, b} is found from the linear fit of logR versus logZ. Equation (27) is an approximation that may be good under a limited set of circumstances. This conclusion becomes obvious by inspection of the DSD moment in Equation (26). Two requirements would need to be satisfied for Equation (27) to be exact: would need to be a constant for all time and drop terminal velocity ( ) would need to be accurately represented by a power law. Neither of these requirements are satisfied, except under very limited circumstances. Nevertheless, Equation (27) can be very useful as long as it is understood that it is just an approximation. Figure 7 shows the Z-R curves generated with a parametric plot by running through a range of values (approximately 1 through 7). Three curves are shown from the FCM-DSD model with three values of , = 8000 (mm −1 m −3 ), = 1, and = 0. The limits of integration D1(w) and D2 as shown in Figure 5, are D1(0) = 0.3 (mm) and D2 = 5.5 (mm). In addition, the National Weather Service (NWS, United States of America) standard Z-R (see [35]) is shown for comparison. Figure 8 shows two sets of curves for comparison with widely varying values of . From these figures it is obvious that the Z-R line shifts with values of , corresponding to a different value of A. There are also accompanying changes in the slope. Changing and also shift the Z-R lines, producing different values of A and b. Setting = 0 turns off the effect has on the FCM-DSD (see Equation (22)). The short dotted and long dotted lines will converge to the solid line for each shown in Figures 7 and 8. This parameter has the effect, on the DSD model of Equation (22), of mixing in a percentage of the standard exponential DSD from 0 to 100%. Since may change during a rainfall event, as well as , a Monte Carlo approach is used. In Figure 9, the FSM-DSD is used to find a {A, b, w} set. Each line is composed of 1000 uniform random distribution Z-R points with ranging from 2000 to 32,000 and ranging from 2.2 to 5, with = 1. The red line is a linear fit to the = 3 m s −1 updraft case, the gray line is the = 0 no vertical wind case, and green is the = −4 m s −1 downdraft. Table 2 summarizes the power law fits. Figure  10 is a plot of Table 2 with the NWS standard Z = 300 R 1.4 and the A-b pair plotted as a yellow square. Since N 0 may change during a rainfall event, as well as Λ, a Monte Carlo approach is used. In Figure 9, the FSM-DSD is used to find a {A, b, w} set. Each line is composed of 1000 uniform random distribution Z-R points with N 0 ranging from 2000 to 32,000 and Λ ranging from 2.2 to 5, with β = 1. The red line is a linear fit to the w = 3 m s −1 updraft case, the gray line is the w = 0 no vertical wind case, and green is the w = −4 m s −1 downdraft. Table 2 summarizes the power law fits. Figure 10 is a plot of Table 2   Verifying the prediction of Figure 10 can be performed by plotting {A,b, } of real data. The A and b can be derived directly from disdrometer data. Vertical wind data can be measured by a collocated microwave wind profiler or by an acoustic SODAR system, if available. Other sources include radiosonde balloon data and global or regional forecast models.  Figure 11 is a plot of the average vertical air motion from z = 0 (or surface) to z = height of 10,000 m  Verifying the prediction of Figure 10 can be performed by plotting {A,b, } of real data. The A and b can be derived directly from disdrometer data. Vertical wind data can be measured by a collocated microwave wind profiler or by an acoustic SODAR system, if available. Other sources include radiosonde balloon data and global or regional forecast models.  Figure 11 is a plot of the average vertical air motion from z = 0 (or surface) to z = height Verifying the prediction of Figure 10 can be performed by plotting {A,b, w} of real data. The A and b can be derived directly from disdrometer data. Vertical wind data can be measured by a collocated microwave wind profiler or by an acoustic SODAR system, if available. Other sources include radiosonde balloon data and global or regional forecast models. Figure 11 is a plot of the average vertical air motion w A from z = 0 (or surface) to z = z re f height of 10,000 m

ECMWF versus Radiosonde
Atmosphere 2019, 10, 77 15 of 21 for different values of the balloon volume V 0 . This curve shows a dramatic increase to the left as V 0 →m T / ρ 0 , the neutral buoyancy point, since any measured vertical motion would then be entirely attributed to an updraft. To the right, as V 0 increases, the ascent rate should increase, such that a nominally measured ascent velocity would be attributed to an increasing downdraft. Reasonable values of w A are on either side of V 0 ≈ 1.3 m 3 .
nominally measured ascent velocity would be attributed to an increasing downdraft. Reasonable values of are on either side of ≈ 1.3 m 3 . Initially, > , implying that the balloon starts out in the turbulent flow regime. As the balloon rises, decreases until it is less than . The point in time where = defines the time t = . In general, the value of will continue to decrease as the balloon continues to rise until the value is = . This is the time defined as when the flow turns to fully laminar. In Figure 11, the curves for and overlap which implies that jumps from to instantaneously. This also results in an instantaneous jump in and . Further analysis has shown that this artificial gap is a consequence of the model and can be mitigated by model modifications, such as reducing the derivative of ( ). For the purposes of this work, these discontinuities in the flow transition region can be mostly ignored. However, future radiosonde analysis models will certainly address this issue.  Table 3 is the result of taking the equivalent of Figure 11 for all cases considered and finding the best fit such that is equal to the ECMWF vertical velocity . Instead of , the equivalent is shown in Table 3. Note that coincidentally, is numerically equal to for = (6/π) 1/2 = 1.38 (m). The values of found this way are constrained to a range from 1.18 to 1.45 (m).  Initially, Re > Re 2 , implying that the balloon starts out in the turbulent flow regime. As the balloon rises, Re decreases until it is less than Re 2 . The point in time where Re = Re 2 defines the time t = t T . In general, the value of Re will continue to decrease as the balloon continues to rise until the value is Re = Re 1 . This is the time defined as t L when the flow turns to fully laminar. In Figure 11, the curves for t T and t L overlap which implies that Re jumps from Re 2 to Re 1 instantaneously. This also results in an instantaneous jump in C D and w c . Further analysis has shown that this artificial Re gap is a consequence of the model and can be mitigated by model modifications, such as reducing the derivative of C D (Re). For the purposes of this work, these discontinuities in the flow transition region can be mostly ignored. However, future radiosonde analysis models will certainly address this issue. Table 3 is the result of taking the equivalent of Figure 11 for all cases considered and finding the best fit V 0 such that w A is equal to the ECMWF vertical velocity w E . Instead of V 0 , the equivalent D b0 is shown in Table 3. Note that coincidentally, V 0 is numerically equal to D b0 for D b0 = (6/π) 1/2 = 1.38 (m). The values of D b0 found this way are constrained to a range from 1.18 to 1.45 (m).   Table 2. This is compared to the A-b coefficients derived from the disdrometer Z-R computations. Finally, the {A, b, w} points from Table 2 are plotted in Figure 12. The yellow squares in this figure correspond to coefficients, A = 300 and b = 1.4 which are the standard values used in the NWS [35]. The circles with black borders are from the Monte Carlo procedure described in the previous section, as shown in Figure 9. Figures 9 and 10 were generated by a uniform distribution of DSD model parameters Λ and N 0 , from Equation (22), with 2.2 ≤ Λ ≤ 5 and 2000 ≤ N 0 ≤ 32,000. In order to line up the w = 0 simulation based on Equation (22) with the NWS squares in Figure 12, the Monte Carlo limits were modified to, with 1.7 ≤ Λ ≤ 6 and 790 ≤ N 0 ≤ 11,000. This has the effect of shifting the A-b points for the specific w cases in the A-b plane.

ECMWF and Disdrometer
1° × 1° ECMWF averaged data to a reference height of 10,000 m, was calculated for the 11 entries of Table 2. This is compared to the A-b coefficients derived from the disdrometer Z-R computations. Finally, the {A, b, } points from Table 2 are plotted in Figure 12. The yellow squares in this figure correspond to coefficients, A = 300 and b = 1.4 which are the standard values used in the NWS [35]. The circles with black borders are from the Monte Carlo procedure described in the previous section, as shown in Figure 9.

Discussion
The goal of this work is to try to deduce information about the character of the vertical air motion from the disdrometer data. This has been motivated by the postulation that the vertical component

Discussion
The goal of this work is to try to deduce information about the character of the vertical air motion w from the disdrometer data. This has been motivated by the postulation that the vertical component of the wind depends on the fall speed of hydrometeors (see [36]). This has mostly been limited to comparing the disdrometer derived Z-R coefficients A and b with any measurement or prediction of vertical wind w. A major difficulty is that other mechanisms may dominate the DSD such as uncorrelated variations of N 0 . Nevertheless, the results of Figure 12 show a remarkable correlation. For the Athalassa disdrometer to ECMWF data comparison, there is a clear segregation of A-b points for w < 0 and w > 0. In addition, the separation line for both A and b matches the NWS point (yellow square). For the model of Equation (22), the correlation with the sign of w is also displayed. The w dependence of the model b parameter dependence is the opposite of the Athalassa data. In the case of the model data, the w = 0 point was forced to line up with the NWS standard values of A = 300 and b = 1.4 point. Based on this relation, it can be concluded that the coefficients and the power series of the Z-R relation are classified according to the vertical velocity, which is one of the main results of this study.
Radiosonde deduced vertical air motion has many potential error sources. A primary error source is initial balloon volume V 0 (or diameter D b0 ). It is essential to tightly control the initial volume or at least be able to measure it before launch. Another source of error is the balloon total mass m T . Again, it is essential to tightly control the mass at launch time or be able to accurately measure the mass. A challenge in measuring inflated mass is the fact that it must not be corrupted by the buoyancy force. Therefore, measuring the deflated mass and the mass of the gas pumped into the balloon from a K-bottle for example, where the weight of the K-bottle is accurately monitored to determine mass of the gas lost from the bottle to the balloon. Even though many other error sources associated with the balloon and radiosonde payload exist, V 0 and m T are particularly sensitive because of the difference term in Equation (1). Small errors in numbers that are very nearly equal can have significant consequences. Other error sources may actually be larger in magnitude, but their effect would more likely be a scaling error.

Conclusions
The FCM-DSD is based on the exponential DSD and is modified by the relationship between drop terminal velocity v(D) and vertical air speed w. The FCM-DSD has a similar appearance to the popular gamma DSD for w < 0 (see [21]). When points {A, b, w} are plotted, as shown by Figure 12, a clear segregation is seen in the A-w plane for both data and model. The points from data are also clearly segregated in the b-w plane, but the model data are on opposite sides of the w = 0 line. Another surprise is that the segregation lines intersecting at w = 0, yield the NWS standard A-b parameters, 300 and 1.4, respectively.
Generalizing any methodology that makes use of radiosonde data of different time periods or diverse localities should take into account issues related to radiosonde data homogeneity [37].
Radiosonde data are routinely assimilated in numerical weather prediction (NWP) models which are subsequently used in the atmospheric data reanalyses. However, NWP models are generally non-hydrostatic, considering vertical velocity as a prognostic variable of the model. However, no specific assimilation on this variable is normally performed [38]. The model derives its own vertical velocity field from the other meteorological fields in its first time-steps of integration. The reason behind is that no observing capability is able to produce some vertical velocity observations which are comparable to the model vertical velocity (at the scale of its mesh). A drastic increase on spatial resolutions of high-resolution NWP models is needed before these models can resolve the clouds and produce some vertical motion which can be compared to, for example, Doppler radar vertical velocity observations [36,39].
Radiosonde data constitute a unique source of information about the vertical atmospheric profile, hence they are essential for the assimilation efforts in generating climate databases, such as the ECMWF reanalyses [28] (and other reanalyses too, e.g., MERRA-2 [40] and JRA-55 [41]). However, vertical velocities are not directly measured from radiosonde observations and, in this respect, directly measured vertical velocities from radiosondes are not assimilated in reanalyses but the vertical velocity field is derived from other observable fields (refer, for example, to the kinematic method used in ERA-Interim). It is worth noting here that interrelationships between reanalysis data and radiosonde timeseries may also be used to homogenize the latter [42] In this work, it was demonstrated that w can be extracted from radiosonde data if initial balloon volume D b0 is accurately measured, along with an accurate measurement of m T . To accomplish this, vertical velocities from radiosonde data were equated to ECMWF vertical velocity fields, thus solving for D b0 with m T set to a constant 0.9 (kg). The resulting values of D b0 are shown to be within a reasonable range of values.
Future work will focus on using laser disdrometer and radiosonde data to further investigate the relationship between disdrometer derived A-b coefficients and vertical air speed w. The disagreement between model and data of the b parameter location in the w-b plane will be investigated in greater detail, but for now it remains inexplicable. Initial balloon volume at the surface before release (m 3 ) v(D) Drop terminal velocity w Vertical velocity (in m s −1 ) w A Average vertical air motion from z = 0 (or surface) to z = z re f height of 10,000 m w c Calculated balloon still air ascent velocity w E ECMWF vertical velocity w r Balloon ascent rate

Appendix B
First three lines of radiosonde data file of 05:00 UTC 24 October 2012.