SURFACE REFERENCE TECHNIQUE REQUIREDFOR A PRIMITIVE RETRIEVAL ALGORITHM OF GPM/DPR

5A.3

Shinta Seto*1), Toshio Iguchi2), and Taikan Oki1)
1) Institute of Industrial Science, University of Tokyo, Tokyo, Japan
2) National Institute of Information and Communications Technology, Tokyo, Japan
1.

INTRODUCTION

The core satellite of the Global Precipitation Measurement (GPM) mission will be launched in2013 with Dual-frequency Precipitation Radar (DPR)
on. DPR is the successor of the Precipitation Radar(PR) on the Tropical Rainfall Measuring Mission(TRMM) satellite. DPR consists of a Ku-band(13.6GHz) radar and a Ka-band (35.5GHz) radar, theformer of which is similar to PR. DPR is expected
to give more accurate estimates in rain drop size
distribution (DSD) and rain rate by taking advantagesof dual-frequency observations. Methods to retrieve DSD and rain rate from dual-frequency observations have been studied for more than 20 years.
Echo power observed by a space-borne precipitation
radar is converted to radar reflectively factor (Zm) at
some ranges and surface backscattering cross section (σ0) at earth surface. The number of observations by DPR is ideally 2n + 2, where n is the numberof ranges. For simplicity, we assume ideally perfectconditions of observation so that many obstacles
such as ground clutter can be ignored.
Meneghini et al. (1992; called ME92) developed a retrieval method for an airborne
dual-frequency precipitation radar. ME92 assumed
that DSD obeys a two-parameter gamma function
and they retrieved 2n DSD parameters from 2n + 2
observations. They applied a surface referencetechnique (SRT) to estimate path integrated attenuation (PIA) by comparing σ0 at a target pixel andσ0 at nearby no-rain pixels (reference pixels). In
SRT, it is generally assumed that the surface conditions between a target pixel and the reference pixels
0
are same in terms of σ , but this assumption is notguaranteed and the difference in surface conditions* Corresponding author address: Shinta Seto,
Institute of Industrial Science, University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo, Japan,
153-8505; e-mail: seto@iis.u-tokyo.ac.jp

often causes biases in PIA. A SRT is applied in thestandard algorithm of PR (Iguchi et al. 2009) andSeto and Iguchi (2007) found that rainfall-induced
0changes in surface conditions affect σ and yieldbiases in PIA and rain rate estimates. Moreover,instantaneous rain rate estimates over land could belargely biased if the land surface type of a target pixel
is drastically different from that of the reference pixels. In case of airborne observations, unstable
attitude of an airborne could change σ0 and cause
bias in PIA and rain rate estimates (Mardiana et al.
2004; called MA04).
MA04 proposed a new retrieval method
without SRT and tested it with air-borne
dual-frequency observations. While PIA is estimated by means of SRT before the retrieval in ME92,
PIA is assumed before the retrieval in MA04. Afterthe retrieval, PIA can be calculated from the retrieved
DSD and is compared with the assumed PIA. Until
calculated and assumed PIAs are judged to be same,the assumed PIA is modified and the retrieval is iterated. The details of methods in ME92 and MA04
will be given in Section 3. MA04’s method showed
a good performance for the air-borne observations intheir paper, but some following researches such as
Rose and Chandrasekar (2005) and Liao and Meneghini (2005) claimed that MA04’s method fails for
relatively heavy rainfall cases.
Adhikari et al.
(2007) also showed cases of failure in MA04’s
method and proposed a modified method by using
differential attenuation (DA) between the two frequencies. As DA can be estimated accurately only
when rain rate is vertically constant or the Mie effectis negligible, the modified method can not be widelyused.
As no methods without SRT are found towork for all cases, the necessity of SRT is strongly
implied. This study examines the reason why
MA04’s method fails for heavier rainfall, confirms thenecessity of SRT, proposes a new retrieval method,and discusses on the accuracy of SRT required toretrieve surface rain rate within a certain error tolerance.
2.

OBSERVATIONS BY DPR

2.1. Assumptions
DPR observes precipitation in nadir direction (with the incident angle of 0 degree) to obtain Zm
at n ranges and σ0 at the surface (Fig. 1). The
depth of each range is constant to be L [km].
Throughout this study, L is set as 0.25 km. Numbers 1 to n are assigned to ranges from the top to thebottom, then the height of precipitation can be given
as n×L. Rain drops are equally distributed in eachrange both vertically and horizontally, and the DSD
follows a gamma distribution function as shown in Eq.
(1).

N ( D ) = N 0 D μ exp[− (3.67 + μ ) D / D0 ],(1)
where D [mm] is drop size and N [mm-1 m-3] isnumber density. N0 [mm-1-μ m-3], D0 [mm], and μ [-]are DSD parameters. In this study, μ is treated as aknown parameter, and N0 and D0 are parameters tobe estimated. Other particles than rain drops canbe totally ignored.
2.2. Theory
Rain rate (denoted by R [mm h-1]) can becalculated from the DSD parameters as Eq. (2).

R = CR ∫∞D =0v( D) D 3 N ( D)dD,

(2)

where CR is a constant for unit conversion (=0.6π×-3-110 ) and v(D) [m s ] is the falling velocity of rain
drop and is dependent on D. According to Gunn
and Kinzer (1949), we employ the following Eq. (3) tocalculate v(D).

v ( D ) = 4.854 × D × exp( −0.195 D ) .

(3)

Real (non-attenuated) radar reflectively
factor (denoted by Ze [mm6 m-3]) and attenuation
coefficient (denoted by k [dB km-1]) are given as Eqs.
(4) and (5).

Ze = CZ ∫∞D =0

k = Ck ∫∞D =0

σ b ( D ) N ( D ) dD

σ t ( D) N ( D)dD

Fig. 1A schematic figure of assumed DPR observation.

where σb(D) [mm ] and σt(D) [mm ] are backscattering cross section and total extinction cross section
and the both are calculated from Mie theory. Ck isequal to 0.01×log10(e), and CZ is given in Eq. (6).
2

λ4 nw + 2
π 5 nw 2 − 1
2

CZ =

2

2

,

(6)

where λ [cm] is the wavelength of microwave, nw [-]is the refractivity index of water, which depends on
the physical temperature of water (denoted by T [K]).
Measured (attenuated) radar reflectivity
factor (denoted by Zm [mm6 m-3]) at the distance of r[km] from the top of precipitation is calculated as Eq.
(7).
r
Z m (r ) = Z e (r ) exp⎛⎜ − 0.2 × ln10 × ∫ k ( s)ds ⎞⎟s=0⎝⎠ , (7)

where s is a dummy parameter of r. Equation (7)can be rewritten in dB unit as Eq. (8).
r

,

(4)

dBZ m (r ) = dBZ e (r ) − 2 ∫ k ( s )ds
s =0

,

(8)

where dBX indicates 10log10(X) for any variable X.
,(5)

Given that r=rs at the earth surface, PIA [dB] can bewritten as Eq. (9).

PIA = 2 ∫

rs

s =0

k ( s )ds

(9)2.3. In discrete form
At a range, dBZe and k are determined
only by parameters of the range, so they can be
written as Eq. (10) and Eq. (11).
dBZe(N0,D0)=10log10N0+F(D0),(10)
k(N0,D0)=N0×G(D0),(11)
where F and G are functions of D0 and the details are
given as below.

{

F ( D0 ) = 10 log10 CZ ∫G ( D0 ) = Ck ∫

∞

D =0

∞

D =0

},(12)σ b ( D) D μ exp[−(3.67 + μ ) D / D0 ]dD

σ t ( D) D μ exp[−(3.67 + μ ) D / D0 ]dD.(13)

On the other hand, Zm at a range can not
be determined by parameters of the range. Whenthe range locates between r=r0-L/2 and r=r0+L/2 (Fig.
2), the average of Zm (denoted by [Zm]) at the range
can be calculated as Eq. (14).
[Z m ] =1
1
Z m (r )dr = ∫Z e (r ) exp⎛⎜ − 0.2 × ln 10 × ∫ k ( s )ds ⎞⎟drs =0⎝⎠
L ∫r0 − L / 2
L r0 − L / 2
r0 + L / 2

r0 + L / 2

r

.(14)

As Ze and k are constant within the range, Eq. (14)can be simplified as Eq. (15).
[ Z m ] = Z m(r0 )[exp (cL/ 2 )- exp (-cL/ 2 )] ,

(15)

cL

where c is equal to 0.2×(ln10)×k. As c and L arepositive, [Zm] is always larger than Zm(r0). When cor L approaches to zero, [Zm] approaches to Zm(r0).
In other words, when the rainfall is weak enough orwhen the width of the range is thin enough, the attenuation from which the range suffers can be regarded to occur along the two-way path between the
top of precipitation and the center of the range.
This attenuation is indicated by dBA as shown in Fig.
2. Generally, the attenuation can be different fromdBA, but it should be larger than dBA[-] and should
be smaller than dBA[+], where dBA[-] (dBA[+]) are
the attenuation along the two-way path between the
top of precipitation and the top (bottom) of the range.
In many previous researches including
ME92 and MA04, the attenuation is set as dBA[+],

Fig. 2Different definitions of attenuation for a
discrete range.

while the attenuation is set as dBA in some researches such as Kozu et al. (1991). No researches are found to set the attenuation as dBA[-].
To compare our results with those of ME92 andMA04, the attenuation is set as dBA[+] also in thisstudy (except for section 6d), therefore dBZm is given
as Eq. (16).
dBZm＝dBZe－dBA[+].
(16)
Please confirm that the following Eq. (17) is hold.
dBA[-]+kL=dBA=dBA[+]－kL(17)
For the discussion later in this paper, we define a
pseudo observation variable dBZe[-2] as Eq. (18).
dBZe[-2]≡dBZe－2×kL=dBZm+dBA[-](18)
Generally, dBZe[α] is defined as Eq. (19).
dBZe[α]≡dBZe＋α×kL(19)dBZe[α] at a range, for any value of α, is independentof the other ranges as dBZe (=dBZe[0]) is.
2.4. Primitive problem
The primitive problem in this study is toretrieve (N0, D0) of n ranges from dBZm of n ranges at
the two frequencies of the DPR. This is a set of 2n
non-linear equations with 2n unknowns. As shownin section 2a, ideal conditions can be assumed.
Moreover, dBZm is assumed to be free from observation errors. SRT can be used to give estimates
PIA, but the accuracy of PIA is not guaranteed.

BM is essentially same with the retrieval method in
ME92.

3.

In iterative backward method (IBM), thesame procedure with BM is applied without SRT. At
first, PIA is assumed arbitrarily, then (N0, D0) aresequentially retrieved from range n to range 1 as BM.
Once (N0, D0) of n ranges are retrieved, PIA is calculated from the retrieved (N0, D0), and the calculated PIA is compared with the assumed PIA. This
comparison is equivalent to check that dBA[-] atrange 1 is zero or not. If the calculated PIA and theassumed PIA are not same (dBA[-] at range 1 is notzero), assumed PIA turns out wrong. By using differently assumed PIA, the same procedure is iterated till the calculated PIA and the assumed PIA
becomes same. In MA04, PIA is assumed to be
zero at the first iteration, then at the second iterationand later, PIA is assumed to be same with the calculated PIA at the previous iteration.

RETRIEVAL METHODS

3.1. Forward method (FM)In forward method (FM), DSD parameters
are sequentially determined from near range (range
1) to far range (range n) without the use of SRT. At
range 1, dBA[-] is assumed to be zero (as no attenuation occurs over the top of precipitation) anddBZm is observed, so according to Eq. (18) dBZe[-2]is calculated at the two frequencies. Here, we havea smaller problem or a lemma to retrieve (N0, D0)from dBZe[-2] at the two frequencies. This lemma isjust a set of 2 non-linear equations with 2 unknowns,
and is indicated by Lemma[-2]. Generally, a lemma
to retrieve (N0, D0) from dBZe[α] of the two frequencies is indicated by Lemma[α].
When theLemma[-2] is solved at range 1, all the variables at
range 1 become known; dBA[+] can be calculated asdBA[-]+2kL. As dBA[-] at range 2 is known to beequal to dBA[+] at range 1, the same process can beapplied to range 2 after range 1. By repeating thisprocess up to range n, (N0, D0) and other variables ofn ranges can be retrieved.

3.2. Backward method (BM)In backward method (BM), DSD parameters are sequentially determined from far range(range n) to near range (range 1) with the use of SRT.
At range n, dBA[+] is no other than PIA, which canbe estimated by SRT, and that dBZm is observed,
then dBZe can be easily calculated according to Eq.
(16). Here, we have a lemma to retrieve (N0, D0)from dBZe at the two frequencies. As dBZe is sameas dBZe[0], the lemma can be indicated by Lemma[0].
When the Lemma[0] is solved at range n, all thevariables at range n become known; dBA[-] can becalculated as dBA[+]-2kL. As dBA[+] at range (n-1)
is known to be equal to dBA[-] at range n, the sameprocess can be applied to range (n-1) after range n.
By repeating this process up to range 1, (N0, D0) andother related variables of n ranges can be obtained.

3.3. Iterative backward method (IBM)3.4. Comparison of methods
As no attenuation occurs over the top of
precipitation, dBA[-] at range 1 should be zero.
Hereafter, it is called upper boundary condition.
The upper boundary condition is satisfied in FM andIBM, but it is not checked in BM. On the other hand,dBA[+] at range n should be equal to the true PIA.
It is called lower boundary condition. The lower
boundary condition is checked neither in FM nor inIBM. From the view point of the boundary conditions on attenuation, FM and IBM are equivalent,
while BM is something different from the other two
methods. According to the categorization of retrieval methods for single-frequency radar (Iguchi etal, 1994), FM and IBM are kinds of initial value
method, while BM is a kind of final value method.
4.

SOLUTIONS OF LEMMAS

Any of the three retrieval methods for DPR
consists of Lemma[α] (α=-2 or α=0). Except forsolving lemmas, no difficulties are involved in themethods. Therefore, in this section, solutions of thelemmas are investigated. As defined, Lemma[α] isFig. 3A function of Fδ(D0) in case of T=300 and

μ=3.
the problem to retrieve (N0, D0) from dBZe[α] at thetwo frequencies. dBZe[α] can be rewritten as Eq.
(20) from Eqs. (10), (11), and (19).
dBZe[α]i=10log10(N0)+Fi(D0)+α×N0×Gi(D0)×L,(20)
where subscript i indicates the frequency; 1 is for
Ku-band and 2 is for Ka-band. Hereafter, the subscript i can be attached to variables and functions
dependent on the frequency.
For α=0, the third term of the right handside in Eq. (20) is zero, and Eq. (20) can be simplified as follows.
dBZe[0]1=10log10(N0)+F1(D0), (21-1)
dBZe[0]2=10log10(N0)+F2(D0). (21-2)
By taking the difference of the two equations in Eq.
(21), the first term of the right hand side is cancelled
as Eq. (22).
dBZe[0]δ=Fδ(D0),(22)
where subscript δ indicates the difference of Ka band
and Ku band, or Xδ=X2-X1 is hold for any variable X.
Eq. (22) has only one unknown parameter D0. Figure 3 draws the function of Fδ(D0) with μ=3.0 and
T=300[K]. Fδ(D0) takes a local maximum when D0
is 0.97mm. As far as investigated in our study andalso in previous studies, Fδ(D0) takes only one local
maximum in case of rain drops. D0 which gives the
local maximum of Fδ(D0) is indicated by D0s. If
dBZe[0]δ is negative, D0 is uniquely determined from
Eq. (22). However, if dBZe[0]δ is positive, two different D0 can satisfy Eq. (22) except for the case of

Fig. 4

Contours of dBZe[-2]2 in case of T=300
and μ=3. Gray region indicates no solutions are given. Dotted line indicatesthat ∂dBZe[-2]2/∂D0=0.

dBZe[0]δ=Fδ(D0s). One of the two D0 is larger thanD0s and another is smaller than D0s, the former ofwhich yields higher rain rates than the latter. In
order not to miss heavy rainfall rather weak rainfall, it
is a reasonable way to assume that D0 is larger thanor equal to D0s. With this assumption (called
D0s-assumption), D0 is always uniquely determined.
When D0 is determined, by substituting D0 into Eq.
(21-1) or (21-2), N0 can be calculated. As shown
above, Lemma[0] has only one solution underD0s-assumption. Still, when dBZe[0]δ is given to belarger than Fδ(D0s) because of some errors (such as
biases in PIA), no D0 can satisfy Eq. (22). In suchcases, Lemma[0] has no solutions.
Next, Lemma[-2] is considered. dBZe[-2]can be written as follows.
dBZe[-2]1=10log10(N0)+F1(D0)－2×N0×G1(D0)×L, (23-1)
dBZe[-2]2=10log10(N0)+F2(D0)－2×N0×G2(D0)×L. (23-2)
And, the difference in dBZe[-2] between the two frequencies is given as Eq. (24).
dBZe[-2]δ=Fδ(D0)－2×N0×Gδ(D0)×L.
(24)
Here, the solution of (N0, D0) in Lemma[-2] can bevisually investigated by Fig. 4, where the horizontal
axis is D0 and the vertical axis is dBZe[-2] . WhendBZe[-2] and D0 are known, N0 can be calculated
from Eq. (24), and dBZe[-2]i (i=1 or 2) can be calculated by substituting N0 and D0 into Eq. (23). In Fig.

4, the contours of dBZe[-2]2 are drawn with μ=3 andT=300K. Gray shade indicates the region where anegative N0 is calculated. On the border line of thegray shade region, N0=0. We can confirm that theborder line is same with the line drawn in Fig. 3.
The solutions of Lemma[-2] correspond tothe crossing points of the line for a given dBZe[-2]δ
(parallel to the horizontal axis) and the contour line
for a given dBZe[-2]2. It is implied by a visual investigation that Lemma[-2] has two solutions and thatone solution is located left to the dotted line, whichcorresponds to ∂ dBZe[-2]2/ ∂ D0=0, and anothersolution is located right to the dotted line. For example, when dBZe[-2]1=40dB and dBZe[-2]2=35dB(dBZe[-2]δ=-5dB), two solutions are found around
D0=1.4mm and D0=2.1mm. In this case, the both
solutions have D0 larger than D0s, therefore two solutions coexist even with D0s-assumption. Generally,
it can be concluded that the number of solutions ofLemma[-2] is two without D0s-assumption and one or
two with D0s-assumption. Exceptionally, Lemma[-2]can have no solutions when biased dBZe[-2]i are
given.
5.

SOLUTIONS OF
CASE STUDIESPRIMITIVE

PROBLEMS:

When
FM
is
appliedwithout
D0s-assumption, a primitive problem has multiple
solutions including the right solution. Even with
D0s-assumption, multiple solutions exist for somecases. When BM is applied with D0s-assumption,
only one solution is derived, but it is not necessarily
the right solution as the PIA estimate may be biased.
In this section, BM is applied with PIA estimates ofdifferent biases in order to evaluate how biases in
PIA affect the estimate of surface rain rate (rain rate
at range n; indicated by Rs). From this analysis, we
can also see the limitations of IBM and FM.
For the above purposes, simple rainfall
cases are tested. In all the cases, DSD parameters(N0, D0, μ) and T are set to be constant over n ranges
for simplicity though it is not physically realistic particularly for T. μ=3 and T=300K and they are known
parameters. (N0, D0) are unknown and the retrieval
methods do not assume that the parameters areconstant over n ranges. PIA estimates have biases
of between -5 dB and 5 dB independently of thefrequencies. D0s-assumption is employed. When
PIA is biased, dBZe[0] may become larger thanFδ(D0s). In such case, any D0 can not satisfy Eq.
(22) and the retrieval has to be terminated. When
DSD parameters are retrieved up to range 1, dBA[-]iat range 1 (hereafter, simply denoted by dBA[-]i) is
calculated. As is explained, the upper boundary
condition that both dBA[-]1 and dBA[-]2 are zero is
generally not satisfied in BM.
In Fig. 5, the horizontal axis is for the biasof PIA at Ku-band (indicated by ΔPIA1) and the vertical axis is for the bias of PIA at Ka-band (indicatedby ΔPIA2). The lower boundary condition is satisfied at the origin (indicated by O; ΔPIA1=0 and
ΔPIA2=0). Solid line indicates that dBA[-]1=0 anddotted line indicates that dBA[-]2=0. The upperboundary condition is satisfied at the crossing pointsof the two lines.
5.1. Necessity of SRT
In case 1, N0=10000, D0=1.5mm, and
n=10. Rain rate is as weak as about 2.1 mm h-1.
Fig. 5a is shown for this case, where it should be
noted that no lines are drawn as the retrieval is terminated in the left upper part. The upper boundary
condition is satisfied only at O, where the lower
boundary condition is satisfied. In this case, theupper boundary condition and the lower boundary
condition are equivalent, and it is implied that IBM
can give the right solution for this case.
In case 2, N0=10000, D0=2.0mm, and
n=10. Rain rate is about 19.1 mm h-1 and is muchheavier than case 1. In Fig. 5b, the dotted line isalmost overlapped with the solid line while the dotted
contour line is also drawn in the right part of thisfigure. It means that the upper boundary conditionis almost or exactly satisfied not only at O. In otherwords, in some cases that PIA is assumed wrongly(the lower boundary condition is not satisfied), theupper boundary condition is almost or exactly satisfied. For such cases, it can not be expected thatIBM gives the right solution.
To investigate the reason of the difference
between case 1 and case 2, supplemental caseswith smaller n are tested. Case 1-1 is same withcase 1, but n=1. As the two lines are overlapped
only at the origin (Fig. 5c), it is suggested that IBM

Fig. 5Figures to check if the upper boundary condition is satisfied or not. On the solid (dotted) line, theupper boundary condition at Ku-band (Ka-band) is satisfied. a) is for case 1, b) is for case 2, c) is forcase 1-1, d) is for case 2-1, e) is for case 2-2, and f) is for case 2-5. In b), part of dotted line isoverlapped with solid line. In d) to f), an open circle indicates the point A.

can give the right solution. When FM is applied forthis case, Lemma[-2] has two solutions; (N0,D0)=(10000, 1.5) and (82082120, 0.626). As thelatter (wrong) solution has D0 smaller than D0s, it is
not satisfied with D0s-assumption. Therefore, FM
as well as IBM has only one solution for case 1-1.
For cases with any n (DSD parameters are same ascase 1) including case 1 itself, FM should have only
one solution. Moreover, considering the equivalence between FM and IBM, it is suggested that IBM
has only one solution for such cases.
Case 2-1 is same with case 2, but n=1.
As shown in Fig. 5d, the upper boundary condition issatisfied at two points: O and another point (indicatedby A). It means that IBM has two solutions (indicated by the solution O and the solution A). When
FM is applied for this case, Lemma[-2] has two solutions; (N0, D0)=(10000, 2.0) and (125862, 1.596).
The former corresponds to the solution O and thelatter corresponds to the solution A.

As D0 of thesolution A is larger than D0s, the solution A is not
rejected by D0s-assumption. As a result, FM as wellas IBM has 2 solutions for case 2-1.
Case 2-2 is same with case 2, but n=2.
In Fig. 5e, the upper boundary condition is satisfied
at two points O and A. The point A in Fig. 5e isexactly same as the point A in Fig. 5d. When FM isapplied for this case, Lemma[-2] at range 1 is sameas of the case 2-1. If the right solution is selected at
range 1, Lemma[-2] at range 2 is also same as of thecase 2-1. The right solution is obtained when theselection is right at the two ranges. If the selection
is right at range 1 but wrong at range 2, the calculated PIA is biased to the point A. If the selection iswrong at range 1, as dBA[+] at range 1 (=dBA[-] atrange 2) is wrongly estimated, Lemma[-2] at range 2
becomes something different from that at range 1.
In this case, as dBZe[-2]i are biased, this Lemma[-2]at range 2 has no solutions.
It can be said that FM has only two solutions O and A for any n larger than 1 (DSD parameters are same as case 2). Generally, the solution A
can be defined that the right solution is selected at
ranges 1 through (n-1) but the wrong solution is selected at range n (as conceptually shown in Fig. 6).
If the wrong solution is selected at any range up to(n-1), Lemma[-2] at the next range has no solutions.
Accordingly, this primitive problem can not haveother solutions than O and A in case 2 (n=10).
Considering the equivalence of FM and IBM, IBMshould have only two solutions in case 2. Although
it is difficult to see in Fig. 5b if the two contour lines
overlap with each other exactly or not, it is now clearthat the two contour lines cross only at the two points
O and A. Still, along the straight line linking O and A,
the upper boundary condition is almost satisfied. If
any PIA on the straight line OA is assumed in IBM, as
the upper boundary condition is almost satisfied, thebiased PIA can be judged to be right unless numerical calculation is done with very high precision.
Here, we acknowledge that IBM has two
kinds of solutions; they can be called “exact solutions” and “approximate solutions”. The exact solutions satisfy the upper boundary condition exactly
and the approximate solutions satisfy the upperboundary condition with small error which is difficult
to be detected in practical numerical calculation.
The approximate solutions become apparent as nincreases. Fig. 5f is for the case same with case 2
but n=5. The two contour lines are close to eachother compared with case 2-1 (n=1) and case 2-2(n=2), still can be distinguished differently from case
2 (n=10). When multiple solutions exist, IBM cannot select the right solution objectively. In suchcases, PIA estimates by SRT are required to distinguish the right solution from other exact and approximate solutions.
5.2. Required accuracy of PIA
The relative error in surface rain rate derived by BM with biased PIA estimates for case 2 isshown in Fig. 7. The relative error is calculated as[Rs-Rs(O)]/Rs(O), where Rs(O) is surface rain rate of
the right solution. When PIA2 is biased, larger errorin Rs is generated than when PIA1 is biased to theFig. 6

A schematic figure of multiple solutions forprimitive problem. Horizontal direction
corresponds to the bias in attenuation.
Closed circles indicate the right solution
and open circles indicate the wrong solution at each range. The profile of solution
A is shown by thick solid lines. At range
2, there may exist other two solutions asshown by gray circles, but not in case 2.

Fig. 7

The relative error in Rs derived by BM withbiased PIA for case 2. Dotted line shows
the straight line linking the points O and A,
where candidates of solutions are assumed to exist.

same extent.
In IBM, by checking the upper boundary
condition, (ΔPIA1, ΔPIA2) is constrained. It can beassumed in case 2 that candidates of solutions exist
on the straight line linking the points O and A.
Therefore, to select unique solution, PIA estimates atboth frequencies are not necessary. If a PIA estimate at either frequency is given by SRT, by combining this estimate with IBM, unique solution can beobtained. This is a modified IBM. Instead of a PIA
estimate at single frequency, the difference of PIA
estimates at dual frequencies (PIAδ≡PIA2-PIA1) canbe used. As was explained in section 1, PIA estimates are subject to biases caused by the change ofland surface conditions. It is expected that PIAδ has
smaller biases than PIA1 and PIA2 as biases in PIA1
and PIA2 are partly cancelled, though it has not been
supported yet by observations. This modified IBM
is similar to a method proposed in Adhikari et al.
(2007), which used DA instead of PIAδ.

Learned from the results of the previous section, thesolution A is focused. Particularly, (N0, D0) at rangen, Rs, and ΔPIAi of the solution A are examined (Fig.
8). Actually, these variables are independent of n,so they can be obtained just by solving theLemma[-2]. Hereafter, if not explicitly notified, (N0,D0) indicates values at range n. An estimate of thesolution X is indicated by putting (X) to the variable.
For example, D0(A) indicates the D0 at range n of thesolution A.

Approximate solutions with relative errorin Rs of 10% (50%) are called solution B (C). For
case 2, (ΔPIA1, ΔPIA2) is (0.045dB, 0.418dB) for
solution B and (0.138dB, 1.697dB) for solution C.
Required accuracy in PIA at Ku-band is much higher
than that at Ka-band. It seems very difficult to satisfy the required accuracy at Ku-band from the experience of TRMM/PR. Although the authors haveno information on possible accuracy of PIA estimatesat Ka-band by SRT, it is assumed to be same as thatat Ku-band in the rest of this paper. It may be possible to satisfy the required accuracy at Ka-band if
the relative error in Rs of 50% is allowed.

holds that D0(A)<D0(O), Rs(A)>Rs(O), ΔPIA1(A)>0,
and ΔPIA2(A)>0. As long as the inequalities hold,the primitive problem is categorized as type-1. The
upper limit of D0(O) in type-1 is indicated by D0x.
When D0(O)=D0x, the solution A becomes identicalwith the solution O. Here, D0x=2.17mm and thecorresponding Rs(O) is 35.4 mm h-1. Type-1 is further divided into type-1a and type-1b by comparing
D0(A) and D0s. When D0(A) is smaller than D0s, by
D0s-assumption, the right solution O is uniquely selected. In this case, the primitive problem is categorized as type-1a. On the other hand, when D0(A)is larger than D0s, the solution A can not be rejectedby D0s-assumption, and the primitive problem iscategorized as type-1b. D0(O) of type-1a is smaller
than that of type-1b, and the threshold of D0(O) between the two types is indicated as D0w. Here, D0w
is 1.77mm and the corresponding Rs(O) is 7.57 mm-1
h . The case 1 is involved in type-1a and the case
2 is involved in type-1b.
When D0(O) becomes larger than D0x, it
holds that D0(A)>D0(O), Rs(A)<Rs(O), ΔPIA1(A)<0,
and ΔPIA2(A)<0. The upper limit of D0(O) to satisfy
the inequalities is indicated by D0y. When D0(O)approaches to D0y, D0(A) diverges to positive infinity.
D0y is estimated to be around 2.67mm and the corresponding Rs(O) is about 167 mm h-1. When D0(O)6.

SOLUTIONS OF PRIMITIVE
GENERAL RESULTS

PROBLEMS:

In previous section, solutions of primitive
problems are investigated for some cases. In thissection, more general cases are investigated to derive robust results.
6.1.The dependence on D0
Same as the previous section, DSD parameters (N0, D0, μ) and T are set to be constantover n ranges. N0=10000, μ=3, T=300K, and D0 isdifferent by cases (between 0.5mm and 3.5mm).

1) Categorization
Here, the primitive problem can be categorized into several types by D0(O). When D0(O) issmaller than D0s(=0.97mm), the right solution can not
be selected as we assume that D0 should be larger
or equal to D0s. In this case, the primitive problem
is categorized as type-0 and is excluded from thisstudy.
When D0(O) becomes larger than D0s, itFig. 8Categorization of the primitive problem with different D0(O). a) D0, b) surface rain rate, c) PIA at
Ku-band, and d) PIA at Ka-band of the right and wrong solutions are shown.

is between D0x and D0y, the primitive problem iscategorized as type-2.
More two types are defined for D0(O) larger than D0y. When D0(O) exceeds to D0y, D0(A)starts to decrease and it holds that D0(A)>D0(O),
Rs(A)>Rs(O), ΔPIA1(A)>0, and ΔPIA2(A)<0.
In
these cases, the primitive problem is categorized as

type-3. The upper limit of D0(O) in type-3 is indicated by D0z. Here, D0z is 2.99mm and the corresponding Rs(O) is 385 mm h-1. When D0(O) isequal to D0z, the solution A becomes identical withthe solution O. When D0(O) is larger than D0z, theprimitive problem is categorized as type-4. In
type-4, D0(A) is smaller than D0(O) and PIA1(A) issmaller than PIA1(O), but the order of Rs(A) andRs(O) and that of PIA2(A) and PIA2(O) are not fixed.
It might not be very important to consider type 3 andtype 4 because rain echoes at near surface aremasked under extremely heavy rainfall.
2) The role of SRT
The role of SRT is examined for each type.
Approximate solutions B and C are defined as insection 5, but if the relative error in Rs(A) is negative,
approximate solutions with relative error in Rs of-10% (-50%) are called solution B (C). If the relative error in Rs(A) is within 10%, the solutions B and
C are not defined. If the relative error in Rs(A) is
within 50%, the solution C is not defined. In Fig. 8c
and 8d, ΔPIA1 and ΔPIA2 of the solutions B and C
are shown as well as those of the solution A.
While SRT is not necessary for type-1a as
the solution A is rejected by D0s-assumption, SRT isnecessary for the other types. In type-1b, bothΔPIA1(A) and ΔPIA2(A) are positive and thatΔPIA2(A) is nearly ten times as large as ΔPIA1(A).
Therefore, it is much effective to estimate PIA by
SRT at Ka-band rather than at Ku-band for the modified IBM. ΔPIA1(A) and ΔPIA2(A) become zero
when D0(O) approaches to D0s or D0x. It becomes
more difficult but less necessary to distinguish O andA, as the difference in Rs between O and A becomes
small. In type 1b, as ΔPIA1(B) is smaller than 0.1dB
and ΔPIA2(B) is smaller than 0.5dB, it seems to bedifficult to distinguish O from B unless SRT would besignificantly improved compared with the TRMM era.
On the other hand, it can be expected to distinguish
O from C as ΔPIA2(C) is larger than 1dB, while
ΔPIA1(C) is smaller than 0.2dB.
Required accuracy in PIAδ is slightly
higher than PIA2. Still, if biases in PIAi are effectively cancelled by taking the difference, using PIAδleads to better estimates in Rs than using PIA2.
In type 2, both ΔPIA1(A) and ΔPIA2(A) arenegative and that |ΔPIA2| is nearly ten times as large
as |ΔPIA1|. Same as type 1, it is much effective toestimate PIA at Ka-band than at Ku-band. When
D0(O) increases, ΔPIA1(A) and ΔPIA2(A) increase
rapidly, therefore it is not difficult to distinguish O
from A. There is some possibility to distinguish O
from B by means of SRT at Ka-band as ΔPIA2(B) is

Fig. 9Categorization of the primitive problem isshown by different background colors(gray for type-0, light pink for type-1a,
dark pink for type-1b, blue for type-2,
green for type-3, and purple for type-4).
a) Rs(O) and b) ΔPIA2(C) are shown as
contours. μ is set to be 3.

larger than 0.5dB, while ΔPIA1(B) is around 0.1dB.
Moreover, as ΔPIA1(C) is near to 1dB and ΔPIA2(C)is larger than 4dB, it can be expected to distinguish
O from C by SRT at Ku-band as well as by SRT at
Ka-band.

In type 3 and in part of type 4, ΔPIA1(A)and ΔPIA2(A) have different signs and |ΔPIA1(A)| isgenerally larger than |ΔPIA2(A)|. It is more effective
to use PIA at Ku-band rather than at Ka-band. If
PIA at both frequencies are available, using PIAδ
may gives better accuracy.

6.2.The dependence on N0
The same simulation with the previous
subsection 6a is done but for different N0 ranging
from 1,000 (dBN0=30dB) to 1,000,000 (dBN0=60dB)
with the step of 1dB, where dBN0 indicates10log10N0.
In Fig. 9a, different types are shown bydifferent background colors on the plane of (D0,
dBN0) with the contour of Rs(O). Note that D0s is
independent of N0 as it is determined by Eq. (22).
D0w, D0x, D0y, and D0z are dependent on N0, and theytend to become small as N0 becomes large. When
N0 is larger than 54dB, type-1a does not exist.
Furthermore, when N0 is larger than 57dB, type-1b is
not found. Except for these high N0 cases, Rs(O) at
the border of different types are not very sensitive to
N0. The upper limit of Rs(O) in type-1a is 5 to 10
mm h-1, that in type-1b is around 20 to 50 mm h-1,that in type-2 is around 150 to 200 mm h-1, and thatin type-3 is nearly 500 mm h-1. Roughly speaking,
SRT is not necessary when Rs(O) is lower than 5 mm
h-1, on the other hand, SRT is required when Rs(O) ishigher than 10 mm h-1.
The contour of ΔPIA2(C) is shown in Fig.
9b. When D0(O) is close to D0x or D0z, as the solution C does not exist, ΔPIA2(A) is shown instead of
ΔPIA2(C). In type 1b, |ΔPIA2(C)| is generally larger
than 1dB where N0 is smaller than 50dB. In type 2,
except for the case that D0(O) is close to D0x,
|ΔPIA2(C)| is larger than 1dB and exceeds to 4dB
when N0 is smaller.
6.3.The dependence on μ
Up to here, the third DSD parameter μ is
fixed to be 3, but the results may be significantly
changed by different μ. In this subsection, μ is fixed
to be 1, and the same simulation with subsection 6bis done (Fig. 10). D0s for μ=1 is 0.71mm andsmaller than that for μ=3. Other thresholds of D0(O)Fig. 10 Same as Fig. 9, but for μ=1. In a), thethreshold below which MA04’s method
can give the right solution is shown according to RC05.
are also smaller for μ=1 than those for μ=3. However, the upper limit of Rs(O) in each type is nearlysame with the case of μ=3. Moreover, the contourof ΔPIA2(C) for μ=1 shown in Fig. 10b is basicallysimilar to that for μ=3 in Fig. 9b. It means that required accuracy of SRT at Ka-band is not very sensitive to μ.
Rose and Chandrasekar (2005) tested

MA04’s method (a kind of IBM) for many cases with
μ=1, and empirically determined the conditions thatMA04’s method could give the right solution. If D0 issmaller than D0* given in the Eq. (25) [as Eq. (41) of
their paper], the right solution is selected by MA04’s
method.
2⎛b ⎞D 0 = ⎜⎜ a + 0.5 ⎟⎟ ,
Nw ⎠
⎝
*

(25)

where Nw is normalized intercept parameter ofgamma distribution, and can be related with D0, N0,
and μ by Eq. (26).

3.67 4 Γ( μ + 4)
N w = N 0 D0
6 (3.67 + μ ) μ +4 . (26)μTwo variables of Eq. (25) a and b were empiricallydetermined and were dependent on n. The lines
corresponding to Eq. (25) are drawn in Fig. 10a forthe cases of n=6, 9, and 12. As values of a and b
were derived under the conditions that Nw is notlarger than 8000 and D0 is not larger than 2.5mm, thelines in Fig. 10a are drawn within in this limit. It isseen that the three lines are drawn in the region oftype 1b. As MA04’s method tends to select thesolution with the smallest PIA among multiple candidates of solutions, it generally selects the solution
O rather than the solution A for type-1b. As it becomes difficult to distinguish O from approximate
solutions when n increases, it is convincing that thelines shift to weaker Rs(O) when n increases.
Adhikari et al. (2007) also tested MA04’s
method and showed (N0, D0)=(8000, 1.1) as a successful case and (N0, D0)=(8000, 1.5) as a failurecase. In their study, it was set that μ=0 and n=17.
We confirmed that the former case belongs to
type-1a and the latter case belongs to type-1b bydrawing the same figure with Fig. 10a but for μ=0(not shown). The previous studies found the limitation of MA04’s method in the region of type 1b.

Fig. 11

Figures to visually examine the solutionsof Lemma[α]. a) α×N0 when dBZe[α]and D0 are given. Shaded area indicates α×N0>0 and plain area indicates α
×N0<0. Contours are drawn for dB|α×
N0|. b) Contours of dBZe[α]2. Dottedline indicates that ∂ dBZe[α]2/ ∂ D0=0.
α=1 for shaded area and α=-1 for plainarea.

6.4.The dependence on the discrete from of dBZm
As explained in Section 2c, dBZm at a
discrete range is defined in Eq. (16) up to here in thispaper. However in this subsection, the attenuationis set as dBA as shown in Eq. (27), to check if different definitions of dBZm affect the conclusions of
this paper.
dBZm=dBZe－dBA.
(27)

With this definition, the second equality of Eq. (18)does not hold. Instead of that, the following two
equations hold.
dBZm＋dBA[-]=dBZe－1×k×L≡dBZe[-1], (28-1)
dBZm＋dBA[+]=dBZe＋1×k×L≡dBZe[1].
(28-2)In FM, dBZe[-1] is given instead of dBZe[-2], then aLemma[-1] has to be solved at each range. In BM(and IBM), dBZe[1] is given instead of dBZe[0], then aLemma[1] has to be solved at each range.
The characteristics of Lemma[-1] andLemma[1] are explained in a general way as below.
Following Eq. (29) can be derived from Eqs. (10),(11), and (19) for any .
dBZe[α]δ=Fδ(D0)+α×N0×Gδ(D0)×L(29)
When dBZe[α] and D0 are given, α×N0 can becalculated. In Fig. 11a, where the horizontal axis isD0 and the vertical axis is dBZe[α]δ, gray shaded area
indicates that α×N0>0 and plain area indicates that
α×N0<0. Contours are drawn for dB|α×N0|. As
N0 should be positive, solutions of Lemma[α] canexist only in gray shaded area when α>0 and only inplain area when α<0. It can be seen that theboundary of shaded area and plain area (α×N0=0)
is identical with the line in Fig. 3. If α=-2 is applied
for plain area, dBZe[-2]2 can be calculated as shownin Fig. 4. Same as this, α=-1 is applied for plainarea and dBZe[-1]2 is calculated. On the other hand,
α=1 is applied for shaded area and dBZe[1]2 is calculated. In Fig. 11b, dBZe[-1]2 and dBZe[1]2 areshown together.
The dotted line indicates ∂dBZe[α]2/∂D0=0 (α=-1 for plain area and α=1 forgray shaded area). From a visual investigation, thefollowings are strongly suggested. Lemma[-1] has
two solutions, and at least one of them has D0 larger
than D0s. It is not always possible to get uniquesolution even with D0s-assumption. Lemma[1] has
unique solution when dBZe[1] is negative and it has
three solutions in some cases when dBZe[1] ispositive. In the latter case, two of the three solutions have D0 smaller than D0s. Therefore, it is always possible to get unique solution of Lemma[1]with D0s-assumption. In addition to the above, it
should be kept in mind that the lemmas have nosolutions in some cases that biased dBZe[α]i are
given.
Figure 12 is same as Fig. 9 but with thedefinition of Eq. (27). D0(O) at the boundary ofFig. 12 Same as Fig. 9, but the Eq. (27) is used todefine the attenuation at a discrete range
in stead of Eq. (16).
different types is slightly larger in Fig. 12a than in Fig.
9a. Corresponding Rs(O) is also slightly larger in
Fig. 12a than in Fig. 9a. ΔPIA2(C) as shown in Fig.
12b is nearly same with that in Fig. 9b. It can beconcluded that the main results of this study are notsignificantly changed by different definitions of dBZm.
7.

SUMMARY AND CONCLUSIONSA primitive problem for the GPM/DPR,
which derives DSD parameters from observed Zm, isinvestigated in this study. This can be treated as amathematical problem to estimate 2n unknowns from
2n equations. Retrieval methods of this problem
such as FM, BM, and IBM, can be decomposed intolemmas to estimate 2 unknowns from 2 equations.
Lemmas involved in FM generally have two solutions,so FM yields multiple solutions for the primitive
problem. On the other hand, lemmas in BM haveunique solution under D0s-assumption, therefore BM
can yield the right solution if the PIA is accurately
estimated by SRT.
IBM is equivalent to FM as the upperboundary condition is checked but the lower boundary condition is not checked in both methods.
There are two types of solutions in IBM: exact solutions and approximate solutions. Exact solutions
satisfy the upper boundary condition exactly andthey are equivalent to solutions derived in FM.
Approximate solutions almost satisfy the upperboundary condition. When n increases, as thenumerical error is accumulated, approximate solutions become apparent.
The primitive algorithm is categorized into
several types. In type-1a, as the right solution canbe selected by D0s-assumption, IBM can yield theright solution. In type-1b, as the solution A can not
be objectively rejected, IBM can not be expected toalways yield the right solution. However, in MA04’s
method (a kind of IBM), which tends to select thesolution with the smallest PIA, the right solution is
selected when rain rate is relatively weak. The
upper limit of Rs(O) in type-1b is generally less than
50 mm h-1 (though it is dependent on N0 and μ), but
the upper limit of Rs(O) when the right solution canbe derived in MA04’s method is much smaller because of approximate solutions.
SRT is necessary to solve the primitive
problem except for some weak rainfall cases. By
combining with IBM, PIA estimate by SRT at either
frequency gives unique solution. In type 1b andtype 2, required accuracy of PIA estimate by SRT ismuch higher at Ku-band than at Ka-band. So, it ismuch effective to estimate PIA at Ka-band if the accuracy of SRT is same at both frequencies. To
estimate Rs with a relative error of 50% or better, therequired accuracy of PIA estimates at Ka-band isaround 1dB in type 1b and around 4dB in type 2.
The above results are not strongly dependent on N0,
μ, and the definition of dBZm in discrete form.
Though some part of this paper is notmathematically rigorous (such that the number ofsolutions is judged by a visual investigation), theresults are believed to be right for most cases. It isnecessary to note that this study assumes many
ideal conditions (such that no observation errors inZm). To compensate such errors, higher accuracy
should be required for SRT. On the other hand, thisstudy does not assume that the relationship between
ranges. As Rose and Chandrasekar (2006) assumes that DSD parameters are expressed by linear
function of altitude, additional constraint may relax
the requirement for SRT.

Acknowledgements. This study is a part of resultsof a collaborative research “Development of surface
reference technique for dual-frequency precipitation
radar (PI: S. Seto)” between the Japan Aerospace
Exploration Agency (JAXA) and the University ofTokyo.
REFERENCES
Adhikari, N. B., T. Iguchi, S. Seto, and N. Takahashi,
2007: Rain retrieval performance of adual-frequency precipitation radar techniquewith differential-attenuation constraint. IEEE
Trans.
Geosci.
RemoteSens.,
45,
2612-2618.
Gunn, R., and G. G. Kinzer, 1949: The terminal velocity of fall for water droplets in stagnant air.
J. Meteorol., 6, 243-248.
Iguchi, T., T. Kozu, J. Kwiatkowski, R. Meneghini, J.
Awaka, and K. Okamoto, 2009: Uncertainties
in the rain profiling algorithm for the TRMM
precipitation radar. J. Met. Soc. Japan, 87(A),1-30.
Kozu, T., K. Nakamura, R. Meneghini, W. C. Boncyk,
1991: Dual-parameter radar rainfall measuremetn from space: A test result from anaircraft experiment. IEEE Trans. Geosci.
Remote Sens., 29, 690-703.
Liao, L., and R. Meneghini, 2005: A study ofair/space-borne dual-wavelength radar for
estimates of rain profiles. Adv. Atmos. Sci.,
22, 841-851.
Mardiana, R., T. Iguchi, and N. Takahashi, 2004: A
dual-frequency rain profiling method without
the use of a surface reference technique.
IEEE Trans. Geosci. Remote Sens., 42,2214-2225.
Meneghini, R., T. Kozu, H. Kumagai, and W. C.
Boncyk, 1992: A study of rain estimation
methods from space using dual-wavelengthradar measurments at near-nadir incidence
over ocean. J. Atmos. Oceanic Technol., 9,364-382.
Seto, S., and T. Iguchi, 2007: Rainfall-inducedchanges in actual surface backscatteringcross sections and effects on rain-rate estimates by spaceborne precipitation radar. J.
Atmos. Oceanic Technol., 24, 1693-1709.
Rose, C. R., and V. Chandrasekar, 2005: A systems
approach to GPM dual-frequency retrieval.
IEEE Trans. Geosci. Remote Sens., 43,
1816-1826.
Rose, C. R., and V. Chandrasekar, 2006: A GPM
dual-frequency retrieval algorithms: DSD
profile-optimization method. J. Atmos. Oceanic Technol., 23, 1372-1383.