Stability of the dual-frequency radar equations and a new method applied to the GPM’s Dual-frequency Precipitation Radar (DPR) data

A new algorithm is proposed that estimates two parameters of the particle size distribution (PSD) at each range bin from the Global Precipitation Measurement’s (GPM’s) Dual-frequency Precipitation Radar data. The equation that expresses the relationship of the PSD parameters between adjacent range bins is derived. By including the attenuation effect within the bin in the discretized equation, the new algorithm alleviates the double-solution problem when attenuation within the bin is sufficiently large. The stability of the solutions to the equation depends on the value of the mean diameter Dm and its gradient with respect to range in the case of liquid precipitation. If the critical diameter above which the dual-frequency ratio (DFR) of radar reflectivity factors becomes a monotonically increasing function of Dm, the backward processing of the equation provides a stable or moderately diverging solution, unlike the forward processing that often gives unstable solutions. To provide a set of initial conditions without using the surface reference technique (SRT) in the backward processing, an initialization method using the Hitschfeld-Bordan attenuation correction method is proposed and tested. The proposed algorithm may provide a tool for investigating the assumptions used in various algorithms.


I. INTRODUCTION
T HE Dual-fruquency Precipitation Radar (DPR) onboard the core satellite of the Global Precipitation Mission (GPM) consists of a Ku-band radar (KuPR) and a Ka-band radar (KuPR) [1]. It was designed to measure backscatter by precipitation in the same scattering volumes at two different frequencies so that it would be able to provide two independent pieces of information about the scattering characteristics of the particles at each scattering volume. The effective radar reflectivities at the Ku and Ka bands differ slightly owing to the fact that the non-Rayleigh scattering effect, or size effect, is more pronounced at the Ka band than at the Ku band. More specifically, if we denote the radar reflectivity factor in dB units at the Ku band by dBZ e1 and that at the Ka band by dBZ e2 , this difference appears in the differential frequency ratio (DFR) that is defined as and Ka bands are also written as dBZ e (Ku) and dBZ e (Ka), respectively in later sections to increase readability.) If we characterize the distribution of precipitation particle size by two parameters, one with the mean diameter D m and the other with a quantity N 0 that is proportional to the number density N t , all radar parameters can be expressed as functions of D m and N 0 provided that other physical and scattering properties of particles are known. Under such conditions, since Z e1 and Z e2 are linearly proportional to N 0 , DFR becomes a function of D m alone. Therefore, once DFR is obtained, D m is expected to be estimated by using the inverse function of DFR(D m ). This logic is the basic idea of estimating D m by a matchedbeam dual frequency radar. Its feasibility for snow parameter estimation was shown by airborne radar data [2]. Nevertheless, after the launch of the GPM satellite, it turned out that estimating D m independently at each range bin was a formidable task. The reasons for the difficulty are multifold. The major difficulty lies in the accurate attenuation correction without the accumulation of biases at either frequency. In order to calculate DFR, we need to estimate Z e1 and Z e2 from measured radar reflectivity factors Z m1 and Z m2 by correcting the attenuations. To correct for the attenuation, we need to know the specific attenuations k 1 and k 2 at each range, but to know them, we need D m and N 0 . Although it is theoretically possible to estimate the attenuation from the nearest range bin to the next successively by solving the equations that relate Z m (r), Z e (r), and k(r) to D m (r) and N 0 (r), this selfreferencing nature of equations tend to accumulate errors in one direction and often gives unreasonable estimates of D m and N 0 . There are several error sources that cannot be totally eliminated. They include noise and fluctuation in measured signals, calibration errors, differences between the assumed particle size distribution (PSD) models and the actual PSD, and other assumptions such as environmental conditions [3].
Attenuation correction by using the path-integrated attenuation (PIA) estimated by the surface reference technique (SRT) [4] is possible and enables an estimation of the DFR at each range bin starting from the last gate of the column and progressing backwards. However, because of other concerns and the difficulty of estimating D m reliably in many cases, the initial operational algorithms (V04, V05 and V06) use a conservative approach in which D m is not estimated independently at each range [5]. The operational algorithms are those used to produce the operational GPM products disseminated by NASA and JAXA. They are called standard algorithms and their products are called standard products in this paper. This paper examines the general properties of PSD parameter estimation and proposes a new algorithm that estimates D m and N 0 at each range. Section II reviews a differential form of the PSD estimation equation and its properties. Section III introduces discretized equations that include the effect of attenuation within the range bin in which the PSD parameters are estimated. By including the attenuation effect within the bin, the double-solution problem is partly mitigated. Section IV describes the general characteristics of the algorithm that include how the dual-solution problem is mitigated. Section V shows a few examples of PSD parameters estimated from DPR data with this algorithm. Section VI discusses the issues of the PSD estimation and Section VII summarizes the paper.

A. Formulation of the problem
We assume rain is uniform in the lateral directions perpendicular to the range direction of radar so that the measurement of rain echo can be treated as a one-dimensional problem. Under this assumption, the expected value of the measured radar reflectivity factor Z m , in dB units (dBZ m ), at range r can be expressed in terms of the effective radar reflectivity factor Z e , in dB units (dBZ e ), and an attenuation factor given as the integral of the specific attenuation k: dBZ mλ (r) = dBZ eλ (r) − 2 r r0 k λ (s) ds, (λ = 1, 2) (2) where k is expressed in dB per unit length, and r 0 is the range at which the radar echo starts. Since we assume that we measure the same rain profiles at two frequencies, any parameters at these frequencies are distinguished by subscript λ where λ = 1, 2. Both Z e and k are functions of the PSD N (D). If N (D) can be well approximated by a model function with two parameters θ 1 and θ 2 , both Z e and k become functions of these two parameters. They can be expressed mathematically as follows: and k λ (θ 1 , θ 2 ) = c k σ tλ (D)N m (D; θ 1 , θ 2 ) dD.
Here, σ bλ is the backscattering cross section and σ tλ is the extinction cross section. The coefficient c Zλ is defined by c Zλ = λ 4 /(π 5 |K λ | 2 ) with K λ = (m 2 λ − 1)/(m 2 λ + 2) where m λ is the complex refractive index of the particle for electromagnetic waves with wavelength λ. c k is a proportional constant that changes with the unit of k λ . Since we use the unit of k λ as dB per unit length as in equation (2), c k = 1/q where q = 0.1 ln(10) is the conversion factor from dB to neper.
Taking the logarithm of (2) and differentiating it with respect to r, we obtain Or expressing it in terms of θ 1 and θ 2 , we obtain If we have radar reflectivity measurements at two frequency channels λ = 1, 2, then we have two equations of the form (6), and we get a simultaneous pair of differential equations for θ 1 and θ 2 .
The inverse of A is given by where |A| is the determinant of A. The PSD, N (D), can be expressed as the product of the scale factor N t and the normalized density function p n (D): We define M i as the i-th moment of N (D): To characterize the PSD by two parameters, they must be a combination of N 0 that is proportional to the total number density and a parameter that characterizes the distribution of p n (D). In this paper, we choose θ 1 = 10 log 10 (N 0 ), and (12) θ 2 = 10 log 10 (D m ). (13) where both N 0 and D m are defined individually by a combination of the third and fourth moments of N (D) as D m is the volume weighted mean diameter, and N 0 is proportional to N t . The intercept parameter N w is related to N 0 by Note that in the case of ice particles, the D m derived from the equations in the paper refers to the unmelted particle diameter whereas in the operational data D m refers to the melted equivalent.
It is known that the natural variations of Z e and k can generally be approximated well by functions of D m for each wavelength once they are normalized by N 0 [6] [7]. In other words, (17) and Then the matrix A and its inverse become where Once a model PSD function is selected, we can calculate f λ (θ 2 ), g λ (θ 2 ) and their derivatives. Since dBZ mλ (r) are measured at discrete ranges r = r i (i = 1, ..., n), they can be used to solve the pair of differential equations (7) to obtain PSD parameters θ 1 (r) and θ 2 (r) numerically. Once these parameters are estimated, we can calculate the rainfall rate from them. This is the basic idea of retrieving two PSD parameters and estimating accurate rainfall rate from matched dual-wavelength radar data.
Unfortunately, however, there are a few problems associated with equation (7). They include a possible double solution problem, instability of the solution, and no solution cases depending on the given data and the assumed PSD model.
It is worth noting here that equation (7) contains the measured data dBZ m,λ by their derivatives only. As a result, constant biases in dBZ m,λ are irrelevant as long as the initial conditions are properly given. In other words, calibration of echo power is important only in the process of determining the initial conditions.

B. Stability of the solution
The solutions to differential equation (7) are determined solely by the initial conditions. This fact indicates the importance of the selection of the initial conditions. If the initial conditions deviate from the true conditions, the obtained solutions may differ from the actual profiles of θ 1 and θ 2 . Whether small deviations of the initial conditions from the true values increase with range or decrease to zero can be analyzed by looking at the stability of the system of equations for the difference of two set of solutions to (7) that result from two initial conditions that differ by a small amount. More specifically, let the solution obtained from the initial condition (θ 10 , θ 20 ) be (θ 1 (r), θ 2 (r)), and a near-by solution obtained from the initial condition (θ 10 + ∆θ 10 , θ 20 + ∆θ 20 ) be (θ 1 (r) + ∆θ 1 (r), θ 2 (r) + ∆θ 2 (r)).
For small (∆θ 10 , ∆θ 20 ), the equation for (∆θ 1 (r), ∆θ 2 (r)) can be obtained by taking the first order terms of the Taylor expansion of (7) and be expressed in the matrix form as follows.
d dr After some tedious calculations, the elements of B turn out to be A prime indicates differentiation with respect to θ 2 . The stability of equation (21) is determined by the trace and the determinant of the coefficient matrix B [8]. It is known that the stationary point, which is ∆θ 1 (r) = 0 and ∆θ 2 (r) = 0 in the current case, is a stable and convergent point (sink) if det B > 0 and Tr B < 0. If det B > 0 and Tr B > 0, it is an unstable and divergent point (source). If det B < 0, it is a saddle point and unstable.
In a case of uniform rain in which θ 1 and θ 2 do not depend on range r so that dθ 2 /dr = 0, the trace and determinant of B are We define the dual-frequency ratio of specific attenuation, DFk, by DFk is a function of θ 2 alone similarly to DFR. Note that DFR = f 1 − f 2 and DFk = q(k 1 /k 2 )(g 1 − g 2 ). Since (f 1 − g 2 )k 2 − (f 2 − g 1 )k 1 > 0, the sign of Tr B depends only on the derivative of DFR with respect to θ 2 . The sign of det B is determined by (g 1 − g 2 )/(f 1 − f 2 ). DFR changes its sign from negative to positive at θ 2 = θ 2c where θ 2c corresponds to the critical diameter D mc at which DFR takes its minimum value. Similarly, g 1 − g 2 changes its sign at the diameter D mk at which DFk takes its minimum value. As a result, when D m increases from 0, the sign of det B changes from positive to negative at D mk , and then to positive again at D mc if D mk < D mc . Figure 1 shows how DFR and DFk as functions of D m change with the temperature and the shape factor µ of the Gamma PSD for liquid water particles. The D mc at which DFR takes its minimum value changes with µ. Different temperatures do not change D mc much, but change the depth of the DFR curve. On the other hand, DFk is relatively insensitive to µ, but the depth and location of the minimum depend on temperature. One other quantity of interest is D ms which is the value of D m above which the inverse function of DFR(D m ) is single-valued. Figure 2 shows how D mc , D mk and D ms change with µ for different temperatures.
In the uniform rain case, if we solve equation (7) backwards, the sign of Tr B changes, but the sign of det B remains the same. From these observations, it can be said that the forward  In the interval [min(D mk , D mc ), max(D mk , D mc )], neither solution is stable (see Fig. 3). It is worthwhile to note that changing the parameters θ 1 and θ 2 to any other related pair of parameters does not change the stability because the linear transformation of infinitesimal parameters does not change the trace and determinant of the coefficient matrix after the transformation.
The mathematical theory of stability can be used to describe the behavior of solution for r → ∞. For the radar problem, however, the range of r is finite. As a result, even if the solution were mathematically convergent, the actual solution obtained with arbitrary initial conditions would not necessarily converge over the finite range. This fact also implies that the solutions do not necessarily diverge from the true solution very rapidly in the cases in which the true solution is not a convergent point. The rate of convergence or divergence is rather complicated; examples of this can be found in section IV.
The behavior of solutions near the stationary point can be analyzed by the eigenvalues and the corresponding eigenvectors. With a Gamma PSD model, the two eigenvalues have a very large ratio of about 100 or more at all realistic values of D m and N 0 . This fact implies that the trajectory of the solution to a stationary point is actually nearly parallel to the eigenvector that corresponds to the larger eigenvalue and does not converge to the stationary point in practice. The trajectories of solutions are effectively the same as in the case with the determinant 0 so that the solution approaches a point on the line of stable fixed points of which direction is determined by the eigenvector for the smaller eigenvalue. Since the angle between these two eigenvectors is rather small, if the equation is solved from the initial values that deviate from the stationary solutions, because of the finite range, the solutions need not end at a stationary point, even in stable cases.
In the non-uniform rain case, because of the additional terms in b 12 and b 22 , the signs of both det B and Tr B depend not only on θ 2 but also on θ 1 and dθ 2 /dr. The additional terms are independent of θ 1 , but the first terms in (23) and (25) are proportional to N 0 (= exp(qθ 1 )), and change significantly with θ 1 . As a result, even a very small deviation of dθ 2 /dr from zero drastically changes the signs of det B and Tr B for a relatively small θ 1 . In fact, a gradient of dθ 2 /dr = ±0.1 dB/km is large enough to modify the original domain of convergence in the uniform case. Figure 3 shows the stable domains for both forward and backward processings when dθ 2 /dr = 0, +0.1 and −0.1 dB/km. It can be said that in forward processing, the stable domain shrinks when there is a gradient in θ 2 regardless of its sign. However, the stable domain increases in the backward processing by adding a new stable region in the negative θ 2 domain when dθ 2 /dr < 0 for relatively small θ 2 . It can be seen in (23) and (25) that if dθ 2 /dr and k λ increase by the same factor, the dependence of the signs of det B and Tr B on θ 2 remains the same. In other words, if dθ 2 /dr increases by a factor of a and θ 1 increases by 10 log 10 (a), the whole patterns in Fig. 3 shift vertically along the N w axis by 10 log 10 (a) and their dependence on θ 2 does not change. For example, if dθ 2 /dr changes from 0.1 dB/km to 1 dB/km, the patterns in Fig. 3 (b) and (c) shift upward by 10 dB. Note that the stability does not depend on dθ 1 /dr, i.e., the gradient of N 0 does not affect the stability.
Although the boundaries of the stable domains change by the gradient of θ 2 , the direction of the eigenvector for the larger eigenvalue remains effectively the same, and the absolute value of the smaller eigenvalue is always very close to 0. Because of these properties, the convergence and divergence patterns of ∆θ 1 and ∆θ 2 in the non-uniform case are similar to those in the uniform case. The eigenvalues of the matrix B determine the convergent and divergent rates of ∆θ per unit distance. Areas with a positive eigenvalue correspond to a diverging area, i.e., the area in which the error in the initial conditions increases with the processing distance. If we accept as a practically stable domain a small eigenvalue area in which the error increases only moderately, the practically stable domain increases substantially. For example, the area with the larger eigenvalue less than q corresponds to the domain in which the increase of the error rate is less than 1 dB/km (≡ 26%/km). The area defined by this condition shares most of the practically important part (dBZ e (Ku)> 20 dBZ) of the stable domain in the uniform rain case. Figure 3

(d) and (e)
show such examples. All these properties of the differential equation (7) are applicable to the discrete equations that will be discussed in the remaining sections.
It is easy to understand that the equation (7) has a singular point at D m = D mc because A −1 has f 1 − f 2 in the denominator which makes A −1 diverge at D m = D mc . Since the right-hand side of equation (7) is discontinuous at D m = D mc , there is no continuous solution of θ 2 that crosses the boundary at D m = D mc . The solution diverges when it hits this boundary. To cross this boundary, we need to use a discretized system of equations with a finite range step. The double value problem and the no-solution issue will be discussed for such discretized equations in the next section.

III. DISCRETIZED EQUATION
A new discretized system of equations is proposed in this section. Instead of discretizing the differential equation (7) by approximating it by finite steps in r, we start with equation (2). We derive a set of equations that relate the PSD parameters in a range bin to those in the adjacent bin. Once such equations are found and if a set of initial conditions are given properly, we can calculate the PSD parameters at all bins.
We assume that the measured data are available only at discrete range bins: r = r i (i = 1, . . . , n, n ∈ N ). r 1 is the range to the center of the first range bin in which the precipitation echo starts. We introduce the following notation. and Here A λ,i is the two-way attenuation to range r i at wavelength λ. With these notations, (2) becomes We assume that we can approximate the attenuation between two adjacent bins by the average of the specific attenuations at the centers of these bins multiplied by twice the distance between them, i.e., Then, by taking the difference of (35) at r = r i and r = r i−1 , and rearranging the terms, we obtain

A. Forward equation
Taking the difference of the two equations for λ = 1 and λ = 2 in (38), we obtain Taking the ratio of k 1,i to k 2,i in (38), we obtain, The substitution of (17) into (41) gives From (42), θ 1,i can be expressed as a function of θ 2,i as Similarly, the substition of (18) into (40) gives Equations (43) and (44) together with equations in (39) constitute the equations for θ 1,i and θ 2,i : Substitution of θ 1,i expressed in terms of θ 2,i in (43) into equation (44) gives a single equation for θ 2,i for given C 1,i and C 2,i . Once this equation is solved for θ 2,i , θ 1,i can be calculated by (43). Since C 1,i and C 2,i are given by the measured data and the PSD parameters at r i−1 , the PSD parameters at r i can be obtained by solving (44) and (43). By repeating this procedure for all i from i = 1, we can estimate the PSD parameters at all range bins.

B. Backward equation
By following exactly the same procedure as in the derivation of the forward equation, we can derive the backward equations.
The backward equations for θ 1,i and θ 2,i are with Equations (45) and (46) are identical to equations (43) and (44) except that the sign of ∆r is reversed and that C 1,i and C 2,i are replaced by B 1,i and B 2,i . B 1,i and B 2,i can be obtained from C 1,i and C 2,i by replacing subscript i − 1 by i + 1 and by changing the sign of ∆r.
A simplified flow chart of the backward processsing is shown in Fig. 4.

IV. BASIC PROPERTIES OF THE ALGORITHM
Since dBZ mλ,i is given at all i (i = 1, . . . , n) for both λ = 1, 2, if θ 1,i−1 and θ 2,i−1 are known, we can calculate dBZ eλ,i−1 and k λ,i−1 and hence C 1,i and C 2,i as well. Similarly, if θ 1,i+1 and θ 2,i+1 are known, we can calculate B 1,i and B 2,i . Substitution of (43) into (44) gives the forward equation of θ 2,i for given C 1,i and C 2,i , and substitution of (45) into (46) gives the backward equation of θ 2,i for given B 1,i and B 2,i . In either case, the equations depend only on two parameters C 1,i and C 2,i , or B 1,i and B 2,i which contain the difference of dBZ mλ,i between two adjacent range bins. As a result, the solutions only depend on the initial conditions and the slope of dBZ mλ,i , and do not depend on the absolute values of dBZ mλ,i . This property is common with the differential form of equations mentioned in section II. Both the forward and backward equations are basically derived from equation (37) together with (35) which can be rewritten as Both sides of this equation give the attenuation estimate at the midpoint between range r i−1 and r i , i.e., at r = r i−1 + ∆r/2 = r i − ∆r/2. In other words, this equation gives the condition that the attenuation calculated forward from the bin center of bin i − 1 and that backward from the bin center of bin i should agree at each frequency.

A. Double-value problem
This method does not solve the problem of dual solutions. In fact, the number of the solutions of equation (44) or (46) change with the value of C 1,i and C 2,i or B 1,i and B 2,i . To explain how the number of solutions changes, we define the left-hand side of the forward equation (44) as L F (θ 2 ; C 1 , C 2 ) or simply L F (θ 2 ) and the left-hand side of the backward equation (46) as L B (θ 2 ; B 1 , B 2 ) or simply L B (θ 2 ) after substituting θ 1 of (43) and (45) into θ 1 in (44) and (46), respectively. When we do not need to distinguish L F (θ 2 ) and L B (θ 2 ), we simply write L(θ 2 ) for brevity. The subscript i is omitted below because the equations hold for an arbitrary range bin. Because of the smallness of the factor DFk(θ 2 ) in (43) and (45), L F (θ 2 ) (L B (θ 2 )) depends much more on C 1 (B 1 ) than on C 2 (B 2 ). In other words, L F (θ 2 ) (L B (θ 2 )) depends much more on the Ku-band than the Ka-band data. Figure 5 shows how L F (θ 2 ; C 1 , C 2 ) (L B (θ 2 ; B 1 , B 2 )) changes with various values of dBZ e (Ku) when true θ 2 = −1 (dBmm) and ∆r = 1 km. C 1 , C 2 , B 1 and B 2 are calculated from dBZ e (Ku) and the given conditions. Note that C 1 and B 1 are approximately equal to dBZ e (Ku). As can be seen from the results, the domain of double solutions of θ 2 increases with dBZ e (Ku) for L F (θ 2 ), but decreases with dBZ e (Ku) for L B (θ 2 ). When dBZ e (Ku) is sufficiently large, L B (θ 2 ) becomes a monotonically increasing function of θ 2 so that the double solution problem disappears. This property can be explained as follows. Equation (45) gives θ 1 as a function of θ 2 . Its dependence on θ 2 is predominantly determined by (45) and hence the contribution of the positive slope from the second term in L B (θ 2 ) overwhelms the negative slope of DFR(θ 2 ) for small θ 2 .
On the other hand, the domain of dual-solution increases with increasing C 1 in the forward processing. Figure 6 shows how the number of solutions depends on the combination of C 1 and C 2 or B 1 and B 2 when ∆r = 1 km. Note that the value of C 1 (or B 1 ) is close to Z e (Ku) because both dBZ m1,i − dBZ m1,i−1 (or dBZ m1,i − dBZ m1,i+1 ) and the attenuation within the bin are generally small. The number of solutions shown in Fig. 6 is the solution within the interval of [−3, 6] for θ 2 . The no solution region in the lower right corner of the graph comes from the upper limit of θ 2 = 6. The right-hand boundary between single and double solutions in the forward processing changes with the lower limit of the allowable domain of θ 2 .
The graph of L B (θ 2 ) shows that the equation L B (θ 2 ; B 1 , B 2 ) = B 1 − B 2 may have up to three solutions because L B (θ 2 ) goes to −∞ as θ 2 goes to −∞. However, in order to avoid unrealistic solutions, the searching domain of the solutions is limited within the interval [θ 2,l , θ 2,u ].
In the examples shown later in this paper, we choose the lower limit θ 2,l = −2 and the upper limit θ 2,u = 6. These numbers correspond to D m = 0.63 mm and D m = 4.0 mm. Within this searching domain, the number of solutions is 0, 1, or 2 depending on the values of C 1 and C 2 in the forward processing and B 1 , and B 2 in the backward processing.
Let us define L d as the minimum value of L(θ 2 ) for θ 2 ∈ [θ 2,l , θ 2,u ] and θ 2,d as the value of θ 2 that makes L(θ 2,d ) = L d (Fig. 7). We also define L s as the value of L(θ 2 ) above which the inverse function has a single solution, and θ 2,s as the corresponding diameter: L(θ 2,s ) = L s . Similarly we define L l = L(θ 2,l ) and L u = L(θ 2,u ). L s = L l but θ 2,l = θ 2,s in general. However, note that θ 2,l = θ 2,d = θ 2,s and L l = L d = L s , when B 1 is sufficiently large and L(θ 2 ) is a monotonically increasing function over [θ 2,l , θ 2,u ].
We take θ 2,l large enough so that L l < L u even if the attenuation within a bin is very large in the forward processing. Under this condition, we can define four possible cases. Let us define ∆C = C 1 − C 2 and ∆B = B 1 − B 2 . Then these four cases are categorized by the interval in which ∆C or ∆B falls: Note that θ 2,d ≈ 0 and θ 2,s ≈ 1 in the case of µ = 3 correspond to values of D m of 1 mm, and 1.26 mm, respectively, and are fairly typical values of D m for light to moderate rain rates.
When ∆B is between L d and L s , two solutions are possible; one in domain I, the left-hand side solution, and the other in domain II, the right-hand side solution. Both solutions satisfy the equation, but they give different estimates of Z eλ and k λ for the following adjacent range bin.
When (44)). In this case,  the obvious choice of the solution is θ 2,d . However, θ 2,d does not satisfy the equation. The difference L B (θ 2,d ) − ∆B (or L F (θ 2,d ) − ∆C) must be accounted for the calculation for the next range bin. Without this correction, the residual error accumulates and the retrieved estimates of the PSD parameters no longer reproduce the measured profiles of Z mλ (λ = 1, 2). In practice, the difference L B (θ 2,d ) − ∆B is subtracted from B 2 in the next step. This compensation method effectively attributes the difference to the Ka-band measured reflectivity factor, Z m2 , and ignores the weak dependence of L B (θ 2 ) on B 2 . The same correction is applied when there is no solution when ∆B > L u .
In the ice region, the DFR is always positive and approaches monotonically and asymptotically 0 dB as D m decreases. In this case, θ 2,l = θ 2,d = θ 2,s and domains I and II disappear. If the measured data gives ∆B < L l ≈ 0, the best choice without any constraint would be θ 2 = −∞ or D m = 0. However, a small value of D m or a large negative value of θ 2 (= dBD m ) gives a huge estimate of θ 1 (= dBN 0 ). This large estimate of θ 1 often causes an overflow in computation of k in which θ 1 appears in the exponent. To avoid this practical problem, we need to define a reasonable minimum allowable D m or θ 2,l . The need for setting a minimum allowable θ 2,l is also applicable to the liquid precipitation.
In contrast to the backward processing, the dual-solution domain of θ 2 increases as the attenuation within the range bin increases in the forward processing. However, since the forward processing is applied only in the solid precipitation region in which the attenuation can be ignored, and since DFR is a monotonic function of D m for solid precipitation, the dualsolution problem does not occur.
The double-value problem is mitigated with the backward processing when the attenuation between adjacent bins is large. To make L B (θ 2 ) a monotonically increasing function of θ 2 , the attenuation term ∆rk 2 in equation (46) needs to be at least one dB. For the DPR, dBZ m,λ is available at range intervals of 0.125 km. With this sampling interval of ∆r, k 2 must be larger than a few dB/km. Such a large attenuation cannot be expected in light rain. In actual situations large attenuations appear only in rather heavy rain. Such heavy rain is generally associated with large D m at which DFR is positive so that the mitigation of the double-solution issue by the proposed algorithm is not effective as long as small ∆r is used. Increasing ∆r will increase the attenuation term and hence the single solution cases. The feasibility of this approach is discussed in the discussion section.

B. Stability
To examine the stability of the solution for the discretized equations (44) and (46), we assume a uniform rain field in which the true Z eλ,i and k λ,i are constant, or equivalently θ 1,i and θ 2,i are constant. Equation (44) for a uniform case (θ 1 = θ 1,i = θ 1,i−1 and θ 2 = θ 2,i = θ 2,i−1 ) becomes where ∆DFRm = DFRm i − DFRm i−1 . Thus, the stationary solution lies on the line defined by the following equation.
The true solution is on this line, but it may not be reached from any initial point in finite steps. Figure 8 shows the flow of the solutions in one step with ∆r = 0.125 km from various initial points in the D m -N w space in the backward processing. The base of each arrow is the starting point, and its tip indicates the solution in the following bin. The directions of arrows are nearly parallel to the line defined by θ 1 = dBZ e (Ku)−f 1 (θ 2 ). Their magnitude decreases as the initial point approaches the stationary solution line defined by (50). Since these two lines cross each other with a small angle, convergence to (or divergence from) the true solution is not apparent. The solutions remain close to the stationary line without approaching the true solution. A closer look reveals that in the unstable case ( Fig. 8(a)), arrows are directed away from the stationary solution line, while in the stable case ( Fig. 8(b)), they are directed toward the stationary solution line as long as the initial points are not far from the true stationary solution. If the initial point is far from the stationary line, the solution after a final step may end up at a point separated substantially from the true solution. As mentioned in section II, since the ratio of the absolute values of the two eigenvalues is very large, a stable solution tends to converge to a point on the stationary line which is not necessarily close to the true solution. These characteristics imply that the selection of the initial conditions is very important. For example, in the case shown in Fig. 8 Nevertheless, if the appropriate initial conditions are properly selected so that they do not cause any significant biases in the estimates of dBZ e,λ , all such solutions satisfy the attenuation condition given by equation (37). This attenuation difference is mainly determined by the Ka-band attenuation. Since the attenuation in the Ka band is highly correlated with the rainfall rate irrespective of the PSD, the estimate of R is generally close to the true value even if N 0 and D m estimates are not close to the true values in the convergent region as long as the total path attenuation is properly given.
The solution from the forward method often gives an unrealistic profile in ranges where the attenuation becomes significant. Nevertheless, when the total path attenuation remains small, the forward method gives rather reasonable profiles. This is the case for snow profiling above the bright band in stratiform rain.
In a solid precipitation region in which the attenuation due to both scattering and absorption is small, we can assume that the attenuations are negligible. In such a case, equation (44) degenerates to DFR(θ 2,i ) = DFRm i . Then, the D m estimates are determined mostly by measured reflectivity factors that are nearly identical to the true effective reflectivity factors. The solutions of PSD parameters in such a region by the forward method effectively depends only on the assumed PSD and scattering models.
This fact can be verified by comparing D m estimated by the forward method with the corresponding D m directly calculated from DFRm. Figure 9 (a) and (b) show the estimated D m at 7 km above the sea level from the forward method and from DFRm measured from the DPR data obtained from orbit 31343, scans 2650∼2749 in the vicinity of hurricane Dorian on 04 September 2019. The agreement is generally very good. Figure 9 (c) shows the averages of D m profiles from the Dorian case. The estimated average D m profiles from DFRm and the forward processing method are effectively identical. Such an agreement is almost always the case in ice precipitation regions above the bright band in stratiform rain.
Even in ice precipitation, discrepancies appear when the measured radar reflectivity factor dBZ m is significant. In such a case, even though the absorption by ice particles is negligible, the attenuation due to scattering causes the difference between dBZ m and dBZ e because of the existence of large particles. Since the PSD of such particles may differ from the assumed PSD, the D m estimates may have large biases. In a heavy convective storm, there may be wet hail stones that further complicate the problem.

C. Initial conditions
The convergence property of the algorithm mentioned in the previous sub-section implies the importance of the initial conditions. In the case of forward processing, we can safely assume that the attenuation to the first range gate is negligible and regard Z m as Z e at both frequencies at the storm top provided that the radar is well calibrated. In this case, we can equate DFRm and DFR from which D m can be estimated there. N 0 can be calculated from Z m (≈ Z e ) and D m provided that we know the phase state and density of the particles in the case of solid precipitation. In the case of liquid precipitation, we may encounter the double solution issue.
In the case of backward processing, we need to specify B 1,i and B 2,i at the farthest range. To do so is equivalent to specifying the pair D m and N 0 from which we can calculate Z e1 , Z e2 , k 1 and k 2 , and hence B 1,i and B 2,i . Instead of D m and N 0 , we can also use the attenuation estimates to the bottom range, since Z e1 , Z e2 can be obtained from Z m1 , Z m2 if the attenuations are known. Once Z e1 , Z e2 are given, we can calculate D m and N 0 and then k 1 and k 2 from them. In particular, if the path-integrated attenuations (PIAs) to the surface at the two frequencies are known, we can use them to specify the initial values of B 1,i and B 2,i . Appendix A gives a formula that specifies the initial conditions by the PIA to the surface when there is some unobservable range between the clutter-free bottom and the surface under the condition that Z e 's at both frequencies are constant there.
In addition to the SRT that gives the PIA estimates, however, a different method is examined to specify the initial conditions in order to investigate possible biases that may exist in the SRT's PIA estimates. In particular, when rain is light, the relative errors in the PIA estimates by the SRT become rather large because their absolute values may be comparable to the fluctuation of surface echoes.
One promising method that is similar to the one used in the operational algorithm is as follows [5]. First, assume a power law k 1 -Z e1 relation for Ku band as k 1 = αZ β e1 where β is a fixed constant. Use the Hitschfeld-Bordan method [9][10] to estimate the attenuation corrected Z e1 (r i ; α) profile of Ku band that depends on the adjustable constant α. Since the assumed k 1 -Z e1 relation defines the relation between D m (≡ θ 2 ) and N 0 (≡ θ 1 ) as described in the first box in Fig. 10, we can calculate the profiles of D m (r i ; α) and N 0 (r i ; α) from Z e1 (r i ; α). Using these D m (r i ; α) and N 0 (r i ; α), Z e2 (r i ; α) and k 2 (r i ; α) can be obtained for Ka band from which we can calculate the profile of the attenuated Z m2 (r i ; α). The final step in the procedure is to adjust α to minimize the difference between the estimated Z m2 (r i ; α) and the measured Z m2 (r i ) near the bottom of the profile. More specifically, minimize the following squared sum: where m is the number of bins from the bottom of the profile over which the difference is minimized. To reduce the error caused by the fluctuation of measured data, it is better to use a large m, but for the estimation of the attenuation at the bottom it is better to use a small m. In the examples shown in this paper, m = 5 is used to estimate the attenuations. Once the best α is determined, use D m (r n ; α) and N 0 (r n ; α) as the initial values at the bottom in the backward processing. We call this method for giving the initial conditions the dualfrequency HB initialization (DHBI) method. A flow diagram of this method is depicted in Fig. 10.

V. EXAMPLES
In this section, retrieved profiles of D m , N w and R from the DPR data using the proposed algorithm are shown. The PSD is assumed to follow the Gamma distribution. Functions f 1 (θ 2 ), f 2 (θ 2 ), g 1 (θ 2 ) and g 2 (θ 2 ) are calculated for a collection of spherical particles with shape factors µ = 0, 1, 2, 3, 4 and water temperature, T , of 0, 10, 20 and 30 • C. To avoid errors associated with measured dBZ m near the noise level, only the profiles that contain dBZ m (Ka) ≥ 20 dBZ over more than 10 consecutive range bins (≡ 1.25 km) are processed [11]. Figure 11(a) shows an example of the vertical profiles of D m retrieved with the proposed backward algorithm with µ = 3 at T = 10 • C in a stratiform area. The D m profile from the operational V06X product is also shown by a dotted line. Three profiles are shown: one with dBZ m smoothed by a nine-point Hanning filter before applying the retrieval algorithm, and the other two without any smoothing to the measured Z m data. The latter two profiles are the left-hand side solution and the right-hand side solution. It can be seen that the fluctuation of the retrieved D m is significant without smoothing. Similar but more pronounced fluctuations appear in the N w and R profiles (Fig. 11(b)). To minimize the fluctuations, the examples shown in the rest of this paper are processed with the smoothed data. The example in Fig. 11 is a relatively rare case in which the left and right-hand solutions differ in multiple bins. However, the overall profiles from these solutions are approximately the same. In this example, liquid phase of precipitation is assumed at all height. Therefore, D m shown in the figure from the current algorithm in and above the bright band may not represent the true diameter. The large increase of D m at about 5 km above the sea level corresponds to the bright band. Nevertheless, it is not very much different from the true unmelted diameter in such a region because the DFR as a function of unmelted diameter is only weakly dependent on the density of ice and water particles. The major error comes from the attenuation correction in such regions that affect the estimate of DFR. Since D m from V06X product gives an equivalent melted diameter, it appears to be smaller than the D m estimates from the current algorithm in solid and mixed precipitation regions.
Note that when the left-hand side solutions of θ 2 are chosen, the corresponding estimates of R become unrealistically large in many cases. Accepting a small value of θ 2 tends to give large estimates of R. This is the reason for defining the lower limit of θ 2,l = −2 (≡ D m = 0.63 mm). When ∆B is between L d and L s , there are two possible solutions. Because the solution on the left-hand side is unstable as proved previously and because it often gives discontinuous profiles as shown later, the examples shown in this paper are the data estimated by choosing the right-hand side solutions all the time unless otherwise mentioned. Figure 12(a) shows the estimated D m at 3 km above the sea level from DPR measurements over hurricane Dorian with the initial conditions selected by the DHBI method proposed in this paper. Figure 12(b) shows the same estimates but with the initial conditions given by the SRT. D m estimates from the operational V06X product are shown in Fig. 12(c) for comparison. Figure 13 shows the estimated R at 3 km above the sea level from the same data sets. The R distributions estimated from the current algorithm with the DHBI method and from the operational V06 product show quite similar distributions of R, although for former estimates are somewhat smaller than the latter.
The comparison between the estimates from the current algorithm and those from the operational product can be seen in Fig. 14. R and D m shown in this scatter plots are right-hand side solutions obtained with the DHBI method in backward processing. The path-integrated attenuations (PIAs) to the surface from the test data are on average somewhat smaller than the PIA estimates from the surface reference technique in this example.
Since the solutions from the proposed algorithm depends on the assumed PSD, their dependence is tested with different µ values of the Gamma PSD. Figure 15 shows how the averaged vertical profiles of the retrieved parameters change with the assumed µ value for pixels with a bright band. Liquid phase is assumed at all heights. The retrieved profiles of dBD m , dBN w , R, dBZ eλ and dBZ mλ (λ = 1, 2) that are within 100 scans (4900 beams) in the vicinity of Hurricane Dorian, and that are classified as having a detectable bright-band, are averaged. The figure shows cases with µ = 0, 1, 2, 3 and 4, together with the output from the V06X operational product. It can be seen that the R profiles with different values of µ do not change much, in particular when the DHBI method is used to select the initial conditions. Similar differences are found in convective cases as well, but in the convective cases, the differences between the estimates from the current algorithm and those from the operational product are larger. In the Dorian case, R estimates with the DHBI method is smaller than the standard product substantially. Figure 15(c) shows average profiles of the PSD parameters from a widespread storm over land. In this case, the current algorithm with the DHBI method still tends to give smaller estimates of R than the standard method. Although the true cause of this difference in R estimates is not clear, one of the major causes can be the non-uniform beam filling (NUBF) corrections that are applied in the operational algorithm, but not in the current algorithm. Nevertheless, the difference in the South American case is not as large as in the Dorian case. This kind of difference between over ocean and over land Fig. 13. Same as Fig. 12, but for the estimated rainfall rate R in dB mm/h. may be universal, and might arise either from different PSD characteristics over ocean and land, or from a possible bias in the PIA estimates by the SRT over the wet surface under rain [12] [13]. Note, however, that the magnitude of differences between the current estimates and the standard product varies even among several over-ocean storm cases that have been tested. It may be worth taking the statistics of such differences to investigate their possible relationship to the type of storm system and characteristic PSD.
These examples show to what extent the assumed µ value affects the estimates. Interestingly, the dependence on the assumed temperature of liquid precipitation is so small that the differences are not noticeable in averaged profiles like ones shown in Fig. 15 when the temperature is changed from 0 • C to 30 • C (results not shown).
Note that D m from the operational product shown in the figure is the melted mean diameter, whereas the D m from the current algorithm gives the actual diameter in the case of ice or mixed phase particles. A consequence of this choice is the large D m from the current algorithm in and near the bright band. The scattering cross sections of liquid particles are assumed even in the bright band in the data shown. The D m profiles show that mixed phase snow aggregates in bright band act as large liquid particles in the radar echoes. Since the profiles shown are obtained by backward processing, a discrepancy in the assumption of precipitation phase in the bright band does not affect the estimates below the bright band.
Another interesting issue is the difference between the Fig. 14. Scatter plots of R and Dm and Ka-band PIA to surface between the estimates from the proposed method with the backward processing (Test) and the corresponding parameters from the operational product (V06). Data include all storm types from the Dorian data. R and Dm data include data between 2 and 4 km above the sea surface. current solutions and the operational product when the initial conditions are set by the PIA from the SRT. The reason for the discrepancy is not clear but may come from a different use of attenuation estimates. In fact, the current method uses the PIA estimates from the SRT in all cases to define the initial conditions whenever the reliability flag of SRT shows that the SRT's PIAs are reliable or marginally reliable. By contrast, the operational algorithm uses a weighted combination of the SRT/PIA and the Hitschfeld-Bordan estimates of attenuation. The application of a non-uniform beam filling correction in the operational algorithm may also contribute to differences, particularly in convective cases.

VI. DISCUSSION
Several issues in the dual-frequency PSD retrieval algorithm are discussed in this section.

A. Noise
The major noise component in the DPR received power is fading noise. Since the log-averaged power is recorded, the signal fluctuates around the mean with a gaussain distribution in dB scale. Since dBD m and dBZ e are approximately linearly correlated, estimates of D m also fluctuate similarly in dB scale if no smoothing is applied to the measured data. If R is calculated directly from the estimated D m and N 0 and expressed in linear scale, its positive deviations are generally larger than the negative deviations due to the nonlinear transformation of dBR to R, and the average of R is slightly larger than the R calculated from the average of dBR (Jensen's inequality). The final estimates of R, D m or N w depend somewhat on at which stage we take the average to reduce the fluctuation errors. In the examples shown in this paper, the measured dBZ m 's are smoothed with a 9-point Hanning filter (i.e., over ±0.5 km). No other smoothing is applied. The choice of the Hanning filter is arbitrary. It is adopted because the smearing of curved profiles such as those near the peak of the bright band is smaller than the flat moving average. Other filters such as the Savitzky-Golay filter that preserves the signal tendency relatively well can be used as well.
The dual-frequency algorithm assumes that echoes at two frequencies are available at all range bins. This condition is not met all the time. In fact, quite often, because of the sensitivity difference and a significant attenuation at the higher frequency channel, meaningful echoes are available only at a single frequency at many range bins. As mentioned before, only the profiles that give more than 10 continuous range bins with dBZ m (Ka) larger than 20 dBZ are processed in the examples shown in this paper.

B. Dual-solution
We show in section II that small errors in initial conditions either shrink or increase only moderately in the backward processing if D m > D mc and dBZ e (Ku) > 20 dBZ. Since the majority of rain cases measurable with the DPR belong to this domain, it is natural to use the backward processing. Unless otherwise stated, the examples shown in this paper are processed with the algorithm that chooses the solution in domain II (θ 2 > θ 2,d , right-hand side solution) rather than in domain I (θ 2 < θ 2,d , left-hand side solution) whenever there are two possible solutions. Stability is another reason for choosing the solution in domain II. If we choose the lefthand side solutions in backward processing, an interesting phenomenon happens. Figure 16 shows a retrieved profile of dBD m and dBN 0 from a pair of synthetic dBZ m,λ profiles. This synthetic data is created by the condition that dBZ e (Ku)=30 dBZ and dBD m changes linearly with range. The broken lines show the true profiles, and the solid lines show the retrieved profiles. The true dBD m varies from domain II to I as the range decreases. The backward retrieved solution oscillates between just above θ 2,l and just above θ 2,s . This phenomenon happens because when the solution is in domain I, the solution at the next step tends to move toward θ 2,l by increasing ∆B. But a small increase of ∆B often results in the case in which the solution exists only in domain III (θ 2 > θ 2,s ) because ∆B becomes larger than L s . Once the solution is in domain III, the algorithm tries to find the next solution with a smaller ∆B and decreases D m . As soon as ∆B becomes smaller than L s = L B (θ 2,s ), two solutions are possible, and the solution in domain I is chosen by rule. This cycle is repeated. As the figure shows, the reproduced dBZ m profiles agree nearly perfectly with the given measured profiles because the true attenuations are given by the initial conditions. The retrieved dBZ e profiles are nearly identical to the original profiles.
This simple example indicates that the reproducibility of Z m profiles does not guarantee the correctness of the retrieved parameters. Furthermore, the uncertainty in each estimate of D m must be rather large when we think about the fluctuation of measured signals and the small gradient of L B (θ 2 ) in domain I and II.  Fig. 17 (c) shows that by the forward processing. The assumed phase state of the particles in the retrieval algorithm is changed from liquid to solid abruptly at H=5.5 km which is slightly above the center of the bright band.
In this kind of light precipitation cases in a stratiform system in which attenuation is minimal, the DFRm method and the forward processing give nearly identical estimates of D m in the solid precipitation region as demonstrated in Fig. 9. The backward processing deviates from the forward solutions or DFRm method in many cases. The D m estimates by the backward processing often end up at the lower bound (θ 2,l ) of the searching domain in the solid precipitation region as can be seen in Fig. 17 (a). This discrepancy results most likely from the fact that the attenuation estimated at the bright band top by the backward processing is biased due to an inappropriate attenuation model through the bright band. Although the estimates are biased, however, the variation of D m with altitude by the backward processing is similar to the forward solutions. The discrepancy of D m in the solid precipitation region between the forward and backward processing cannot be resolved no matter the height at which the assumed phase is changed from liquid to solid in the backward processing. It seems necessary to introduce a layer in the bright band in which the specific attenuations at Ku and Ka bands take a rather different ratio than that given by a simple liquid particle model. Note that there are many solutions that end up at the maximum diameter of 4 mm at low altitude in the forward processing because of the divergent tendency of the forward solutions in the liquid precipitation region. It is also worth noting that the left-hand side solutions in domain I in the forward processing do not diverge to the minimum allowable value of θ 2,l = −2 (D m = 0.63 mm), unlike the left-hand side solutions in the backward processing.
The estimates of D m in these figures are obtained by assuming the density of snow as ρ = 0.2 g/cm 3 . As mentioned before, these estimates are relatively insensitive to the assumed value of ρ. However, the corresponding N w estimates depend on the value of ρ because Z e for a given D m depends on ρ. Therefore, even though the discontinuity of the estimated unmelted D m at H=5.5 km is not very large, the corresponding N w estimates show a very clear discontinuity at this height (not shown). To estimate continuous and realistic N w , we need further information or assumptions in the phase transition layer and above.
As mentioned in section III, the proposed algorithm does not solve the dual-solution problem entirely. Nevertheless, when the attenuation within the bin is large enough, L B (θ 2 ) becomes a monotonically increasing function of θ 2 and the dual-solution problem disappears. In fact, with the Dorian data shown in the paper, there are about 50 thousand pair of adjacent bins in rain region that satisfy the conditions for testing the algorithm. Of these 50 thousand adjacent pairs, if we use ∆r = 0.125 km, function L B (θ 2 ) becomes a monotonically increasing function in about 9% of the pairs. In the remaining 91%, L B (θ 2 ) is not monotonic, but the actual solutions are in the single solution domain III (θ 2 > θ 2,s ) more than 80% of the cases and the dual-solution (θ 2 < θ 2,s ) occurs in only about 10% of the total cases. Increasing the range bin size will increase the attenuation between adjacent bins and lower the condition for L B (θ 2 ) to be a monotonic function. If we use ∆r = 0.5 km, for example, the number of cases with a monotonically increasing L B (θ 2 ) increases to 30%, and the dual-solution cases decrease to about 7%. If we use every pair of dBZ m separated by 8 range bins so that ∆r = 1 km, the number of cases with a monotonically increasing L B (θ 2 ) increases to 47%. In the remaining 53%, 41% of solutions are in the single solution domain III, no solution case (∆B < L d ) occurs in 7% of solutions, and only 4% of total cases are subject to the dual-solution problem. The latter number decreases to 1% if the threshold of processing Ka data is increased to 25 dBZ from 20 dBZ. In the South American case, these numbers are 58%, 38%, 3%, and 1% so that the double-solution cases become very rare. Figure 18 shows CFADs of D m from the Dorian and South American cases when ∆r = 1 km is used. Because of the increase of the single solution cases, the left-and right-hand side solution cases show nearly the same distributions in both cases. Judging from their distributions in rain regions below 4 km, we can observe that D m is in domain II rather than domain I in the majority of cases when ∆B < L s . In addition to the stability issue, this observation is another reason for selecting the right-hand side solutions as the solutions of this algorithm in the examples in this paper. However, increasing the step size may violate the assumption of rain uniformity and may result in an adverse effect. In fact, ∆r = 1 km is too large to estimate the PSD parameters in and around the bright band where the specific attenuation at Ka band can change significantly over a short distance.

C. Uncertainties in PSD and scattering model
If valid Z m data are available at two frequencies at n range bins, we have at most 2n independent pieces of information. This fact implies that we can estimate only up to 2n parameters. In the current case, we choose D m and N 0 at n range gates. All other parameters that affect the retrieval must be assumed in some way. For example, they must be defined as  functions of D m and N 0 , or given by environmental data independently provided. In the current case, we need to assume the functions f 1 (θ 2 ), f 2 (θ 2 ), g 1 (θ 2 ) and g 2 (θ 2 ). These functions depend on the PSD and scattering properties of individual particles. The latters depend on the phase, melting fraction, density, shape, and temperature of precipitating particles. In particular, since the specific attenuation depends on the phase state substantially, phase identification becomes crucial for accurate estimation.
Liao et al. [2] show that the snow density does not strongly affect the dependence of DFR on D m when the DFR is taken between X and Ka bands, although the DFR changes somewhat with the shape factor µ of the assumed Gamma PSD. These characteristics of DFR dependence on the snow density and the shape factor are the same for the DFR taken between Ku and Ka bands. Therefore, unmelted D m values estimated by assuming ρ = 0.2 g/cm 3 do not change significantly even if we assume a different ρ value. However, the melted D m and N 0 depend on the assumed value of ρ.
As to the estimated D m for rain, it depends slightly on the assumed temperature of drops because the minimum values of both DFR and DFk depend on the temperature as depicted in Fig. 1. However, since the differences due to the temperature difference are small except near the minimum point of those functions, the estimated D m does not change much even if we change the assumed temperature from 0 to 30 degrees. (D m at the minimum of DFR changes from 1.00 to 1.02 mm when the temperature changes from 0 • C to 30 • C if µ = 3.) The estimated D m depends more on the assumed value of the shape factor µ. (D m at the minimum of DFR changes from 0.78 to 1.07 mm when µ changes from 0 to 4 when the temperature is 10 • C.) This dependence on µ is more problematic because the actual PSD changes substantially and the dependences of k and Z e on D m and N 0 deviate from those derived from any PSD model functions. As a result, even if the true rain is vertically uniform and the attenuation correction is properly applied to determine the initial conditions, the estimated PSD profiles may result in some bias in its absolute values and its vertical gradient.
In a region of dry snow in which attenuation can be effectively ignored, however, the issue is not very severe. We can estimate the unmelted diameter D m from DFRm as mentioned before. The property mentioned in the previous paragraph may be used to identify the existence of hail and graupel by looking at retrieved D m and Z e values as long as the particles are dry. However, if some wet hail exists, the problem becomes formidable. Large hail stones may add multiple scattering echoes that make the problem even more difficult. How to deal with such cases is an interesting issue, but outside the scope of this paper.
We have shown that the right-hand side solutions in the backward processing is stable or only moderately divergent as long as D m is larger than the critical diameter of about 1 mm and estimate D m in rain region reasonably well. However, the estimates within the bright band may not be accurate because of uncertainties in the scattering characteristics in the transition layer from solid to liquid phase. To estimate D m within and above the bright band in the backward processing, appropriate PSD and scattering models are required.
Both non-uniform beam filling (NUBF) and multiple scattering change the apparent dependence of k and Z e on D m and N 0 . In particular, the effective k is no longer a local function of D m and N 0 , but depends on the three-dimensional distribution of precipitation. As a result, attenuation correction becomes a formidable task and reliable correction is virtually impossible when these effects are significant. Since k cannot be expressed in terms of D m and N 0 locally, the current formulation would require some modifications in such cases.

D. Initial values
As mentioned in the algorithm section, the solutions of the proposed method depend critically on the initial values. A self-consistent DHBI method to estimate the attenuations to the bottom bin is proposed and tested. This initialization method, however, depends on the assumed vertical profile of phase state, although its effect is not significant because the attenuation from solid precipitation range to the total path attenuation is generally rather small. Nevertheless, there is a room for improvement. There may be a better way to find initial values of D m and N 0 at the bottom in the backward processing. If estimation with a fixed PSD model and other simplifying assumptions such as no NUBF or multiple scattering effect does not give reliable estimates, improvement of the initialization method alone will not further improve the estimates. Although better ways to specify the initial conditions have not been pursued, further work along these lines is warranted.

E. Other possibilities
In this paper, we discussed a method for estimating θ 1 and θ 2 as functions of range by solving equation (2) with (3) and (4) in terms of differential equations or discretized difference equations. In these methods, solutions are derived from one range to the next successively. There is a possibility, however, that 2n unknowns θ 1,i and θ 2,i can be determined by minimizing the difference beween measured reflectivities and calculated ones without assuming any continuity: where dBZ mc,λ denotes the calculated measured radar reflectivity factor. This method does not work in practice because of the following reasons. First, there are many local minimum points in the 2n-dimensional space of θ 1,i and θ 2,i . The double-value problem associated with DFR also contributes to this problem. It is not easy to find the global minimum point by any numerical methods even if the model assumptions are perfect. Even when there is no noise, the global minimum point is not a sharply defined minimum in some directions in the 2n-dimensional space. Existence of noise in actual data further complicates the distribution of minimum points and may often result in profiles of θ 1 and θ 2 that are rather different from the true profiles. Additional continuity constraints that allow only small changes of the parameters between adjacent bins reduce the acceptable minimum points for the solution, but do not eliminate the multiplicity of solutions. These methods also share the same issues mentioned in this section.
VII. CONCLUSION General properties of retrieving two PSD parameters from matched beam dual-frequency radar data are reviewed and a specific PSD retrieval algorithm for the DPR is proposed. An equation is derived that expresses the relationship of the PSD parameters between the adjacent range bins including the attenuation effect within the bin. This algorithm estimates the mean diameter of particles by solving the equation in one step from the measured radar reflectivity factors at two frequencies and the PSD parameters in the adjacent bin. In the past attempts for estimating the PSD parameters, the attenuation correction and the estimation of D m are carried out in different steps so that a recursive processing is needed in order to include the attenuation effect within the bin itself.
The stability of the solutions to the equation is examined. It is shown that the stability depends on the location of the true D m in the case of liquid precipitation. A small gradient of D m with respect to range may alter the stable domains of solutions substantially from the constant D m case. Nevertheless, if D m is larger than the critical diameter D mc , which is about 1 mm, the backward processing provides either a stable solution in which the initial errors decreases, or a solution whose error increases only marginally with the processing distance. In contrast to the backward prosessing, the forward processing gives stable or moderately diverging solutions for D m < D mc , and diverging solutions for D m > D mc . As a result, errors in the forward processing tend to increase substantially with range in majority of rain profiles measurable with the DPR.
It is also shown that when the attenuation within the bin is sufficiently large, the backward function L B (θ 2 ) becomes a monotonically increasing function of θ 2 so that the doublevalue problem disappears. However, if we use the range step size of 0.125 km, the double-value problem remains in many cases in the DPR data. Therefore, the double-value problem cannot totally be avoided. If the echo profiles are relatively uniform, by increasing the step size, the double-value problem is mitigated in most cases for the signal levels detectable by the DPR.
As a method to provide a set of initial conditions in the backward processing, an attenuation estimation method based on the HB attenuation correction method is proposed and tested. This approach gives reasonable initial conditions and the average of estimated R agrees reasonably well with that estimated with the PIA from the SRT. Since the selection of the initial values has a decisive impact on the solutions, further study may be needed to find a better initialization method.
It is also shown that if the attenuations to the bottom bin are well estimated, the retrieved profiles of R are nearly independent of the shape factor of the Gamma PSD model, even though the estimated D m and N 0 profiles change with the µ value. If the NUBF correction is excluded, the current algorithm produces reasonable R estimates that mostly agree with the corresponding estimates from the operational products on average, even though there are some differences in the estimates of D m and N w . We find some bias in the case of Hurricane Dorian. Since the conditions adopted in the test runs shown in this paper are different from those in the operational processing due to the rather complicated conditions used in the latter algorithm, some differences are not surprising. Nevertheless, the current algorithm may provide a tool to investigate the validity of various assumptions used in the operational algorithm. It may also be used to study the PSD in other more complicated cases such as convective storms with hail and graupel.

CONDITIONS
In the backward equations, we need to give appropriate initial conditions at the farthest range. If the PIA is available at two frequencies, we can give such initial conditions. We assume that Z eλ (λ = 1 for Ku and λ = 2 for Ka band) is constant from the farthest bin (r n ) to the actual surface (r s ). Let the distance from the center of the clutter-free bottom bin to the actual surface be ∆r s : where PIA λ,s is the path-integrated attenuation to the surface for frequency λ. Since we assume that Z eλ is constant between r n and r s , PIA λ,s − A λ,n = 2∆r s k λ,n = 2∆r s exp(q(θ 1,n + g λ (θ 2,n ))) (A.55) Substitution of (A.55) into (A.54) gives dBZ eλ,n = dBZ mλ,n + PIA λ,s − 2∆r s k λ,n (A.56) or in terms of θ 1,n and θ 2,n , θ 1,n + f λ (θ 2,n ) = dBZ mλ,n + PIA λ,s − 2∆r s exp(q(θ 1,n + g λ (θ 2,n ))) (A.57) Therefore, replacing ∆r i by 2∆r s and B λ by B λ = dBZ mλ,n + PIA λ,s , (λ = 1, 2) (A.58) we can find the solutions at range n.
From 1985 to 2019, he was with the National Institute of Information and Communications Technology, Koganei, Japan. He visited the NASA/Goddard Space Flight Center, Greenbelt, MD, USA, from 1991 to 1994, performing the US-Japan collaborative experiment for measuring rain using airborne radar. He is currently a Visiting Research Scientist within the Department of Earth System Science Interdisciplinary Center, University of Maryland, College Park, MD, USA. His work has been mainly related to remote sensing of environment, in particular, issues related to the remote sensing of precipitation from space. He has served on the science teams for TRMM and GPM since their inception and has developed operational radar algorithms for both radars.
Robert Meneghini Robert Meneghini is an employee at NASA Goddard. With scientists from the Communications Research Laboratory of Japan, he co-led a series of airborne field campaigns from 1985 to 2000, the data from which were used to develop algorithms for the TRMM and GPM spaceborne weather radars. He has served on the science teams for TRMM and GPM since their inception and has developed operational radar algorithms for both radars. His research interests include scattering theory, rain and snow retrievals from active and passive instruments and the use of radar for wind estimation over ocean.