A methodology to analyze precipitation profiles using the Tropical Rainfall Measuring Mission (TRMM) Microwave Imager (TMI) and precipitation radar (PR) is proposed. Rainfall profiles are retrieved from PR measurements, defined as the best-fit solution selected from precalculated profiles by cloud-resolving models (CRMs), under explicitly defined assumptions of drop size distribution (DSD) and ice hydrometeor models. The PR path-integrated attenuation (PIA), where available, is further used to adjust DSD in a manner that is similar to the PR operational algorithm. Combined with the TMI-retrieved nonraining geophysical parameters, the three-dimensional structure of the geophysical parameters is obtained across the satellite-observed domains. Microwave brightness temperatures are then computed for a comparison with TMI observations to examine if the radar-retrieved rainfall is consistent in the radiometric measurement space. The inconsistency in microwave brightness temperatures is reduced by iterating the retrieval procedure with updated assumptions of the DSD and ice-density models. The proposed methodology is expected to refine the a priori rain profile database and error models for use by parametric passive microwave algorithms, aimed at the Global Precipitation Measurement (GPM) mission, as well as a future TRMM algorithms.
Precipitation measurements from space play an important role in observational studies of global climate. A significant advance in this area has been made by the Tropical Rainfall Measuring Mission (TRMM), launched in November 1997. The TRMM satellite carries both active and passive microwave sensors to measure rainfall: the precipitation radar (PR) and TRMM Microwave Imager (TMI). The PR is the first spaceborne radar for observing precipitation, while the TMI has a design that is similar to the Special Sensor Microwave Imager (SSM/I) on board the Defense Meteorological Satellite Project (DMSP) satellite, except for improved spatial resolution and an additional pair of 10.65-GHz channels to improve the capability of heavy rainfall detection. These two sensors, individually, or in combination, allow for a greater variety of approaches to examine rainfall distributions. The initial PR and TMI products, however, had significant rainfall bias relative to each other. Version 5 of the TRMM operational algorithms had a 24% discrepancy in the global mean precipitation between the TMI and PR level-2 rainfall products, denoted 2A12 and 2A25, respectively (Kummerow et al. 2000). In addition, the products also showed regional and temporal disagreements in response to the tropical climate variability associated with large-scale circulation and El Niño–Southern Oscillation (ENSO) (Berg et al. 2002; Robertson et al. 2003). Preliminary indications are that the overall bias is being drastically reduced in version 6 of the products. Regional and temporal biases, on the other hand, have been only slightly reduced. Continuous efforts are under way to improve the operational algorithms with the goal of reconciling any inconsistencies.
A number of investigations have been carried out that examine the key factors that are responsible for the sensor-dependent biases in rainfall estimations. Viltard et al. (2000) conducted a case study to show that a common drop size distribution (DSD) assumption for simulating radar reflectivities and microwave brightness temperatures yields a roughly consistent solution of precipitation profiles. Masunaga et al. (2002) focused on the precipitation water content (PWC) and precipitation water path (PWP) derived from the TMI and PR operational products, and found a PR excess in near-surface PWC over TMI in the midlatitude winter and a TMI excess in PWP over PR in the Tropics. They indicated that these biases, combined with a difference in the conversion from precipitation water to rain rate between the algorithms, account for the known discrepancies in the operational products. Ikai and Nakamura (2003) identified several specific error sources that are inherent in each of the TMI and PR algorithms, such as the freezing-level assumption in the TMI algorithm and the Z–R and k–Z relations used by the PR algorithm. Robertson et al. (2003) suggested that systematic changes in DSD that are not accounted for in the 2A25 algorithm may be responsible for additional biases observed during ENSO events.
The above studies suggest the need for a combined radar and radiometer analysis that is grounded in a common physical basis, instead of simply striving to reconcile the independent operational TRMM products. A methodology of such an analysis is proposed in this study, aimed to improve the a priori databases used by the TMI operational algorithm, as well as parallel research efforts (Mugnai et al. 1993; Kummerow and Giglio 1994; Smith et al. 1994; Bauer et al. 2001; Kummerow et al. 2001). These schemes have historically employed only cloud-resolving models (CRMs) to build a priori databases and are, therefore, not necessarily consistent with observed rainfall structures. These schemes can all be improved by constructing databases closer to the rainfall structures that are actually observed by PR. Shin and Kummerow (2003), who did some earlier work in this area, saw this as a critical step to build parametric algorithms that could be employed uniformly across a number of independent radiometers.
A significant improvement has been made since Shin and Kummerow (2003), although the present study follows some of their technique without major changes (e.g., the interpolation of the nonraining parameters across raining areas). Shin and Kummerow (2003) made the first attempt to demonstrate the performance of a parametric retrieval algorithm based on synthetic data. For this purpose, they simplified the procedure to construct the a priori databases, focusing more on testing the applicability of the databases. This paper, in contrast, is mostly devoted to describing the technical aspect of the database development. The database construction scheme has been significantly refined since Shin and Kummerow (2003), so that physical consistency is ensured in the treatment of radar and radiometer measurements. In particular, Shin and Kummerow (2003) used TMI brightness temperatures only to assess the algorithm performance, while the present study directly incorporates brightness temperatures in the analysis to optimize the solution by adjusting the DSD and ice-density models, as described later.
There are some differences in objectives between this study and other work related to combined radar and radiometer algorithms, such as Olson et al. (1996) and Grecu et al. (2004), and the TRMM combined PR and TMI (2B31) algorithm (Haddad et al. 1997). These efforts were concerned primarily with developing operational algorithms. Our goals are summarized as follows. First, this work is aimed at the improvement of a priori databases for use by parametric retrieval algorithms as an extended study initiated by Shin and Kummerow (2003). Second, the current method is also targeted at enabling refinements of CRM through joint efforts with the model community. CRMs describe the dynamical evolution of precipitating cloud systems under given initial and boundary conditions. The legitimacy of those prescribed conditions, however, can be constrained only by observations. A large number of observations, together with nonraining parameter retrievals, gives us an opportunity to investigate rainfall variability in terms of the large-scale environment characterized by SST, humidity, wind speed, and other environmental parameters. It can be useful to examine the relationship between those geophysical parameters and CRM profiles by tracing the PR-selected profiles back in the original database. We see this as a quantitative measure to test CRM in the context of various large-scale environments, as done by a GCM superparameterization technique. Our long-term purpose is, thus, not only to develop a new scheme for precipitation retrievals, but to seek a way to exploit rainfall observations for improving models as well.
Finally, the algorithm architecture is designed so that the known inconsistencies in the TRMM rainfall retrievals between PR and TMI can be directly addressed. A number of studies mentioned earlier in this section have shown that those inconsistencies are closely related to errors in the assumptions built in the algorithms and their biases associated with climate variability. The analysis scheme presented here preserves the PR retrieval with given assumptions without adjustment separately from that with the optimized assumptions using TMI measurements. This approach enables us to explicitly quantify the impact of TMI on the assumptions in the PR-only algorithm, in contrast to other combined algorithms that are aimed at finding solutions that simultaneously invert the TMI and PR data. The present method first employs only PR measurements. To the extent that PR measurements contain explicit assumptions related to DSD and ice particle properties, these are then adjusted to be consistent with TMI observations as well. Uncertainties in PR-retrieved rainfall, such as the DSD assumption, are explicitly verifiable by examining the inconsistency with TMI measurements. This keeps the error propagation transparent throughout the analysis, and provides the basis for eventually interpreting differences between radar and radiometer products. The representativeness of the a priori database, corresponding to the prior probability distribution in the Bayesian framework, is also critical to ensure unbiased retrievals. Once the a priori database is constructed offline, existing Bayesian inversion schemes may be used to perform operational rainfall retrievals.
The applicability of these databases to parametric algorithms, as outlined by Shin and Kummerow (2003), makes such an approach more important in the future as increasing numbers of microwave radiometers with various designs, for example, the Special Sensor Microwave Imager Sounder (SSMIS) and Advanced Microwave Scanning Radiometer (AMSR), are launched and operated. Although the TRMM is currently the only satellite that carries both a radar and a radiometer, the a priori database that is constructed based on the TRMM analysis can be shared by any radiometers that are unaccompanied with a radar. This is a great advantage for the Global Precipitation Measurement (GPM) mission that is being planned for launch around 2010, where a core satellite carrying a dual-frequency radar and a TMI-like radiometer is supposed to algorithmically calibrate constellation satellites with only radiometers aboard. An a priori database constructed from the core satellite provides such a homogeneous resource that all constellation radiometers can employ in their estimates.
The present paper is a detailed description of the methodology used to construct an a priori database from radar, radiometer, and CRM data (section 2). Section 3 discusses the performance of an iterative procedure used to adjust the radar retrieval assumptions to minimize inconsistencies with radiometric measurements. Brief concluding remarks are provided in section 4.
2. Data and method
This section offers a detailed description of the methodology and datasets that are applied to the analysis. The algorithm flow is outlined in Fig. 1. For it to be reproducible, the theoretical background and technical aspects of the analysis are covered in some detail.
a. Satellite sensors and data
The TRMM PR and TMI data are utilized as observational inputs. The TRMM satellite was launched in November 1997 into a sun-asynchronous orbit with a 35° inclination and an altitude of 350 km. The TMI, a conically scanning sensor with a viewing angle of 52.8° at the earth’s surface, has nine channels, as described in Table 1. The TMI swath is 759 km wide with 104 (208) effective fields of view (EFOVs) per scan for channels 1–7 (8 and 9). The down-track interval between adjacent scans is 13.9 km. The PR is a single-frequency (13.8 GHz) radar having horizontal and vertical resolutions of 4.3 and 0.25 km, respectively, at nadir. The PR has a cross-track scan ranging between ±17°, resulting in a swath width of 215 km. The PR sensitivity or minimum detectable echo is 17 dBZ after the system noise is subtracted (TRMM PR Team 2000). The spatial resolutions and swath widths have changed slightly since the TRMM orbital altitude was raised from 350 to 402 km in August 2001. The present study uses data acquired prior to the orbit boost.
The analysis requires four datasets from the TRMM standard data products: the 1B11 TMI brightness temperatures, 1C21 PR reflectivity profiles, the 2A21 PR two-way path-integrated attenuation (PIA), and the 2A25 PR rainfall products (shown in the upper-left corner of Fig. 1). The 1B11 brightness temperatures at channels 1–7 are interpolated so that these channels have the same footprint allocation as channels 8 and 9. The 1C21 dataset contains radar reflectivities without attenuation correction and the rain flag that defines the presence or absence of rain signals. A given PR pixel is designated as “rain certain” when echoes above the noise level are detected in clutter-free ranges. In this analysis, the 250-m vertical resolution of PR is reduced to a 0.5- to 1-km resolution in order to match CRM simulations needed to supply those parameters that are not directly observed by PR (discussed in section 2b). The 2A21 algorithm adopts the surface reference technique (SRT) to estimate PIA, exploiting a depression in the surface scattering of radar echoes in raining areas compared to surrounding rain-free areas (Meneghini et al. 2000). The 2A25 rainfall product is not required as an input by our algorithm. It is used simply for comparison with the results of the present analysis. The 2A25 algorithm derives rainfall profiles from given radar reflectivities with a hybrid scheme combining the Hitschfeld–Bordan equation (Hitschfeld and Bordan 1954) and SRT for the attenuation correction (Iguchi et al. 2000). Because the SRT provides an independent piece of information in addition to radar reflectivities, the 2A25 algorithm allows the coefficients in the k–Z relation, and, thus, the Z–R relation, to be optimized when PIA information is deemed reliable. Iguchi et al. (2000) call this technique the α-adjustment method. Because the k–Z and Z–R coefficients are related to the DSD assumption, the α-adjustment method is philosophically identical to adjusting the DSD model.
b. CRM database
The rainwater and precipitation ice can be derived directly from PR measurements, albeit with some assumptions. For the database construction described here, however, a complete set of atmospheric parameters is needed. This includes those parameters that are not directly measured by PR, such as the temperature structure, water vapor, and cloud water within the precipitating column. To that end, a precomputed set of CRM profiles is introduced and matched with PR reflectivity and PIA, when applicable. This study adopts the source set of CRM simulations used in the latest version of the Goddard profiling (GPROF) algorithm (Kummerow et al. 2001). The database, constructed from the Goddard Cumulus Ensemble (GCE) model (Tao and Simpson 1993) and the University of Wisconsin (UW) Nonhydrostatic Modeling System (NMS) (Tripoli 1992), currently consists of 32 snapshots, including different evolutionary stages of tropical mesoscale convective systems (MCSs), scattered convection, hurricanes, warm/cold frontal rain systems, and extratropical cyclones. Each snapshot is organized in the format of a 136 × 136 to 256 × 256 horizontal grid, with a horizontal resolution of 1.33–4 km and 28 vertical layers with the resolution of 0.5 km below 10 km in altitude and 1 km above 10 km up to 18 km. Precipitation profiles consist of the following five hydrometeor species: rain, snow, graupel, cloud water, and cloud ice. Atmospheric background parameters are also included in the CRM outputs.
c. Raining parameters
The raining parameters or precipitation liquid/ice water content are determined principally by PR measurements (indicated by yellow-colored panels in Fig. 1). The largest sources of errors for precipitation retrievals for a single-frequency, nonpolarimetric radar, such as PR, are the uncertainties in rain DSD and in the microphysical properties of frozen and melting particles. Our analysis starts with the DSD assumption, which is consistent with the 2A25 DSD model used to derive their Z–R relations (Iguchi et al. 2000). The ice-particle characterization is initially the same as provided by the original CRM database. Under these initial assumptions, the CRM profile that best reproduces a given series of radar reflectivities is selected for the solution. The PIA is also incorporated to optimize the DSD model when a PIA signal is deemed reliable by the 2A21 product. By repeating this procedure over all of the raining near-nadir PR pixels, one obtains the complete spatial structure of the raining parameters in the PR-observed domain. A detailed description of the procedure is provided in the following sections.
1) Raindrop size distribution
The DSD model is prescribed in terms of the median volume diameter D0 as a function of the surface liquid water content. Use of D0 is intended to make our DSD model verifiable by ground measurements because it is a commonly used variable to characterize DSDs in polarimetric radar measurements (Bringi and Chandrasekar 2001). This study assumes DSD to be a gamma distribution (in mm−1 m−3),
where N0 is in mm−1 m−3, D0 in mm, and μ is the nondimensional exponent of the gamma distribution that is assumed to be 3, following the PR 2A25 algorithm (Iguchi et al. 2000). For a gamma DSD, the liquid water content (or the mass density of rainwater; in g m−3) W and the radar reflectivity (in mm6 m−3) Z are expressed in terms of N0, μ, and D0 as
where ρw = 1 g cm−3 is the density of liquid water and Γ is the gamma function. Rain rate (in mm h−1) R is derived similarly as
where the drop fall velocity (in m s−1) υ(D) was given by
For a realistic range of D0, W/R are weakly, but nonlinearly, dependent on D0 and μ. Although R is a more conventional quantity that is used to evaluate the rain intensity, we adopt W instead of R as an independent variable for modeling D0 [see Eq. (11) below], because R suffers from an additional ambiguity due to the uncertainty in υ(D). Eliminating N0 by combining Eqs. (2) and (3), one has
Equation (8) then leads to a unique relationship between D0 and W if Z is given as a function of W. A Z–W relation is introduced as
with which a relationship between D0 and W is obtained by eliminating Z in Eq. (8),
The coefficients aw and bw, derived from the same DSD models as are assumed by the 2A25 algorithm, are shown in Table 1 of Masunaga et al. (2002). We adopt the values that are computed for 0°C rainwater, that is, aw = 0.001 998 06 and bw = 0.613 42 for stratiform rain, and aw = 0.003 917 52 and bw = 0.578 55 for convective rain.
The initial DSD model D0,i in the present study is defined as
where D0,s (D0,c) is calculated by Eq. (10) using stratiform (convective) aw and bw. In Eq. (11), D0,i is designed to stay close to D0,s for light rains, but shifts closer to D0,c as W increases so that it reflects the general trend that light rain is produced mainly by stratiform precipitation, while deep convection is responsible for heaviest rainfall. Figure 2 shows D0,i, D0,s, and D0,c as a function of W, as well as four additional DSD models parallel to D0,i by the intervals of ±0.3 and ±0.6 mm. These additional models are employed to adjust the DSD when PIA is large, and will also be used to update the DSD model to minimize the inconsistency between the measured and computed brightness temperatures as described in section 3.
The backscattering and extinction coefficients are calculated based on the CRM-provided liquid water content, along with the DSD model specified by Eq. (11) for use by the radar echo simulation. Mie calculations are applied neglecting a small departure from the sphericity of raindrops.
2) Frozen and melting particles
Frozen and melting hydrometeors have very complicated and highly variable crystal structures. This makes it extremely difficult to model their microphysical and dielectric properties precisely. The irregularity in shape of ice/melting hydrometeors leads to an additional error source in forward radiative transfer simulations that are used by the retrieval algorithms. Melting particles, in particular, are radiatively important because they account for radar bright bands as well as an increase in the attenuation of radar echoes within a melting layer. Microwave brightness temperatures can also be significantly influenced by the enhanced extinction by melting hydrometeors, but a large uncertainty remains because the physical process of particle melting is not well understood (Bauer et al. 1999; Olson et al. 2001; Battaglia et al. 2003).
The particle size distributions of snow and graupel are assumed to follow an exponential form of
and Wi (in g m−3) is the ice (snow or graupel) water content provided by CRMs. The particle density ρi is 1 × 105 g m−3 for snow and 4 × 105 g m−3 for graupel, and the intercept parameter N0,i is given as 1 × 105 mm−1 m−3 for snow and 4 × 103 mm−1 m−3 for graupel. The present study relies on the Maxwell–Garnett model, assuming ice matrix with air inclusions to calculate the effective dielectric constant of frozen particles (Bauer et al. 1999; Olson et al. 2001). The volume fraction of air inclusions is given as 1 − ρi/ρsi, where ρsi is the solid ice density or 0.917 × 106 g m−3. Mie calculations are used, assuming an effective dielectric constant, to compute the backscattering and extinction coefficients for those particles, as done for raindrops. Melting particles are treated by defining a 500-m melting layer below the 0°C height. This domain is divided into a number of sublayers to take into account different stages of the melting process. The melting particle size distribution is linearly interpolated from the snowflake size distribution to the rain DSD, while the dielectric constant of melting particles is kept constant at a fifty–fifty mixture of snow and rain within the melting layer. Finally, the computed radiative properties are averaged over the sublayers. Although this treatment is highly simplified, the resultant extinction coefficient was found to fall onto a moderate value, compared to various models of melting particles shown by Fig. 4 of Battaglia et al. (2003).
3) Radar echo simulation
For a given set of the radiative properties of hydrometeors and the atmospheric gas absorption coefficients, radar reflectivities and PIA are calculated for every CRM profile stored in the initial database. Each CRM profile is accompanied by 25 separate simulations with five DSD assumptions shown in Fig. 2, and by five different assumptions of the ice-density model as mentioned later in section 3. The appendix shows the equations to be calculated, Eqs. (A7) and (A11), as well as their derivations.
Simulated radar reflectivities and PIA constitute the database to be searched in the matching process, where a vertical sequence of simulated radar reflectivities are compared, one at a time, with a ray of PR reflectivities. As such, the physical consistency that is assured in the CRM framework is preserved with respect to the vertical structure of precipitation.
Radar echo simulations are performed only along a nadir view to avoid the geometrical complexity arising from the scan-angle dependence. This constrains the current analysis to be applicable only to near-nadir PR scans. The PR scans in use are ±10 (±5) scans around nadir, or approximately 90 (50) km in swath width, and are correspondingly 7° (3.5°) in scan angle, for a case study (for daily statistics) shown later in this paper. A deviation from nadir is smaller than 1% in terms of the direction cosine of the PR viewing angle for both the cases.
4) Matching procedure
The procedure to identify the optimal rain profile that matches a given PR measurement is described here. The TRMM SRT (or 2A21) algorithm defines the “reliability factor” as the ratio of PIA to the standard deviation in the normalized scattering cross section of the rain-free surface, 〈σ0 (NR)〉 (Meneghini et al. 2000). The SRT estimation of PIA tends to be unreliable when PIA is small (or rainfall is light) and/or when 〈σ0 (NR)〉 is very noisy in the case that neighboring rain-free areas are poorly sampled. The PIA is incorporated in the matching process when the reliability factor exceeds 3 (hereafter denoted as the heavy-rain case). The PIA is neglected (the light-rain case) otherwise. The best-fit solution is defined as the CRM profile with the simulated radar profile that minimizes
where yobs, ysim, and 𝗦ε are the observed variables, simulated variables, and the measurement error covariance, respectively. In the light-rain case, observations and model profiles are compared in terms of radar reflectivities without PIA, that is,
where Zi is the reflectivity at the ith layer and Nl is the total number of layers. The assumed DSD model is fixed to be Eq. (11) (the thick, solid line in Fig. 2). Because rain rates at different layers are not independent of each other (Coppens et al. 2000), 𝗦ε would be a nondiagonal matrix. In the present study, however, off-diagonal components are ignored for brevity, leaving a future study to fill them out. The ith diagonal component of 𝗦ε, wii, is given by
and i0 denotes the layer that contains the 0°C height. The system noise of PR is about 14 dBZ for PR (TRMM PR Team 2000). The modeling errors are attributed to ambiguity in the dielectric modeling of hydrometeors, which is minimal for raindrops and worse for melting particles, as discussed previously in this section. Currently, the modeling errors for raindrops (σrainmdl), melting particles (σmeltmdl), and ice hydrometeors (σicemdl) are given as 0.3, 3, and 1 dBZ, respectively. Because of the limitation of our present knowledge on the hydrometeor dielectric properties, these numbers are only tentative, which qualitatively take into account a larger uncertainty in the dielectric modeling of ice and melting particles.
The error covariance matrix is assumed to be diagonal, as before. Nonzero components in 𝗦ε are same as in the light-rain case, except for the additional diagonal component corresponding to PIA,
By virtue of an additional degree of freedom given by PIA, the best-fit solution is now sought among the entire set of candidate profiles with variable DSD models (all solid lines in Fig. 2) to constrain the DSD assumption. This approach is conceptually parallel to the α-adjustment method of the 2A25 algorithm (Iguchi et al. 2000), although technical details are different.
d. Nonraining parameters
Atmospheric water vapor (WV), cloud liquid water (CLW), surface wind (SW), and sea surface temperature (SST) in rain-free areas are designated as nonraining parameters in the present study. These parameters are not detectable by PR but are needed for the computation of microwave brightness temperatures later in the analysis (section 2f). The nonraining parameters are derived from TMI brightness temperatures for the footprints that are diagnosed as rain free by collocated PR signals, and are then interpolated across all raining and nonraining scenes at the PR resolution. The computational procedure, outlined by blue-colored panels in Fig. 1, is described in this subsection.
Column water vapor (CWV), liquid water path (LWP), SW, and SST are adopted from the TMI products provided by Remote Sensing Systems (RSSs) in the current analysis. Daily datasets are employed for all of these parameters except SST, for which 3-day-averaged datasets are used. The temporal variation in SST is generally negligible on a time scale of 3 days. The daily products are virtually equivalent to the instantaneous retrievals, except for some of the following minor differences: 1) the distributed products are projected on the 0.25° × 0.25° global grid, and 2) two or more snapshots observed at separate times are combined together where different orbits overlap on the same day. Therefore, the original resolution is recovered by extrapolating the gridded data on TMI footprints, and the latitudinal ranges that are higher than 30° are excluded from the current analysis to avoid the orbital overlaps. The RSS TMI products contain two different versions of surface wind speed depending on the channel used for retrieval (10 or 37 GHz). We adopt the 37-GHz version for surface wind speed because it is less susceptible to the sun glitter than the 10-GHz version.
1) Identification of rain-free TMI footprints
Figure 3 schematically illustrates a TMI footprint with a PR pixel to be collocated. The geodesic direction vector pointing at a PR pixel from the TMI footprint center dPR, and the direction from the subsatellite point to the TMI-footprint center dTMI define an angle θ between them as
where dTMI dPR represents the scalar product of dTMI and dPR. Assuming a TMI footprint to be an ellipse with the diameters equivalent to the cross-track and down-track instantaneous fields of view (IFOVs; Δξ and Δη, respectively), the radius of the TMI footprint along dPR, rfp is given as
A PR pixel is considered to fall in the given TMI footprint when
If none of the PR pixels that satisfy this condition are considered to be rain certain, as defined in the 1C21 dataset, the TMI footprint under consideration is defined as rain free.
We utilize the IFOVs of the 19-GHz channel as criteria common to all of the TMI channels for identifying clear footprints. Use of the 19-GHz IFOVs guarantees the clear condition Eq. (26) to be satisfied at all other channels, except at 10 GHz. The 10-GHz brightness temperatures are generally insensitive to the nonraining parameters, except for SST, which uses the 3-day-averaged value and is, thus, available regardless of the presence of rain.
The development of microwave algorithms to derive the nonraining parameters over the ocean has a long history in the satellite remote sensing community (e.g., Chang and Wilheit 1979; Greenwald et al. 1993; Wentz 1997). The present analysis utilizes the RSS TMI data archive, as mentioned previously. The vertical distributions of WV and CLW are determined from the RSS-provided CWV and LWP. For brevity, WV (in g m−3) is assumed to be exponential as a function of z under a given CWV (in kg m−2),
where the water vapor scale height hWV is fixed at 2.3 km. Cloud layers in rain-free areas are assumed to be
where CLW is in g m−3 and LWP is in kg m−2.
Our computed brightness temperatures from the nonraining parameters (as defined above) show good agreement with TMI measurements, except for a small overestimation (underestimation) of brightness temperatures for small (large) values of CWV. This bias results from the difference in radiative transfer schemes and/or optimization procedures used by the RSS retrieval algorithm. We, therefore, made a slight change in CWV from the RSS-provided CWV, CWVRSS, as
to ensure consistency between the nonraining parameters and the brightness temperatures from which they are computed. This modification ensures agreement between computed and measured brightness temperatures at 10-, 19-, and 21-GHz channels for clear-sky areas, while it results in a slight overestimation of computations at the 37-GHz channel.
e. Interpolation of geophysical parameters over PR pixels
Once the nonraining parameters are retrieved from TMI measurements, they are interpolated across the whole observed scenes at the PR resolution. The surface wind SW, at the nth PR pixel SWn, is evaluated as a weighted average computed from surrounding rain-free TMI footprints,
with a Gaussian weight as a function of the distance between the nth PR pixel and the ith TMI footprint dni,
The variance σ2 represents the typical scale of spatial correlation in SW, but the present analysis assumes an arbitrary value of σ = 30 km because of the lack of our knowledge on the physics governing σ. The SST is interpolated in the same manner as for SW.
In contrast to SW and SST, raining PR pixels yield their own WV and CLW, obtained through the PR-matching technique that is described in section 2c. These are used as well for the interpolation in conjunction with the rain-free TMI footprints. The equivalent of Eq. (30) for WV and CLW becomes
where q(z) denotes the vertical profiles of WV and CLW, and ΔTMI (ΔPR) is the mean area coverage represented by a single TMI footprint (PR pixel), defined as the product of the along-scan pixel interval and the down-track scan interval for TMI (PR) or 4.6 km × 13.6 km (4.3 km × 4.3 km). Equation (32) is applied only to rain-free PR pixels.
The interpolated nonraining parameters, together with the PR-derived raining parameters, provide a complete three-dimensional structure of the geophysical parameters at the PR resolution within the central portion of the PR swath. Microwave radiative transfer calculations are then applied to these parameters to examine their consistency with observed brightness temperatures.
f. Brightness temperature computations
The geophysical parameters that are derived from PR observations are not necessarily guaranteed to accurately reproduce TMI brightness temperatures. We argue two major factors that could cause a disagreement in the brightness temperature space. First, the DSD uncertainty still remains after the radar matching process because heavy rainfall that leads to a robust estimate of PIA is relatively infrequent. The initial estimate of DSD is, therefore, not updated for a majority of precipitation. An error in the DSD assumption would result in a bias in computed brightness temperatures through the PR’s misestimation of the rainwater content. Second, frozen hydrometeors are not well constrained by PR because frozen and melting particles are less weighted in the least squares fitting [minimizing Eqs. (14)], and are, hence, less well captured by PR echoes than raindrops are. Therefore, it is not until TMI brightness temperatures, particularly at 85 GHz, are incorporated in the analysis that ice particles aloft are given a satisfactory observational constraint.
Brightness temperatures are computed by a radiative transfer scheme based on the Eddington approximation (Kummerow 1993) with a given set of the extinction coefficient, single scattering albedo, and asymmetry factor, derived in a consistent manner, with the backscattering and extinction coefficients applied to radar echo simulations (section 2c). The atmospheric temperature is defined assuming the constant lapse rate of 6.5 K km−1 bound by SST at the bottom of the atmosphere. Radiative transfer calculations are performed over every single PR-viewing column under a modified plane-parallel approximation that takes into account the intersection of a slant TMI sight line across several adjacent columns. The slant-intersection effect is largest at the 85-GHz channel because of the smallest footprint size and the sensitivity to ice particles aloft. As illustrated in Fig. 4, TMI could see an appreciable amount of ice aloft on the sight line originating from a rain-free surface (arrow with a solid line in Fig. 4) or, on the contrary, could miss the ice, while the associated rainfall below is detected (arrow with a dashed line). A horizontal displacement due to the slant-intersection effect is roughly estimated as zice tan θs, where zice is the altitude of the scattering ice particles and θs is the sensor viewing angle at the surface, and is estimated as 13 km for zice = 10 km and θs = 52.8°, or the TMI viewing angle. The “modified” plane-parallel approximation is designed to properly deal with the slant-intersection effect to ensure a pixel-by-pixel matching in the brightness temperature space.
Once brightness temperatures are computed at the PR resolution, they are convolved with a two-dimensional Gaussian beam pattern,
where ξ and η denote the cross- and down-track coordinates, respectively, Δξ is the cross-track IFOV (Table 1), Δη is the down-track IFOV, and Tb is the brightness temperature at the TMI resolution.
g. Initial case study
In this final part of section 2, the methodology that is introduced above is summarized and applied to a case study in order to briefly demonstrate the performance of the algorithm. A tropical MCS over the central Pacific Ocean that is observed by TRMM (orbit number 13183; 12 March 2000) is used for the case study.
Input data, required in advance, are a set of radar and radiometer observation datasets, an a priori database constituted of CRM-computed precipitation profiles, and a radar profile database containing reflectivities and PIA simulated from the CRM profiles (orange and light-blue ellipses at the top of Fig. 1). First, the precipitation liquid/ice water content is derived from PR measurements. The simulated profile that best fits a given radar measurement is sought and identified, where the PIA is used as well as the radar reflectivity in the matching process to constrain the DSD assumption when robust PIA signals are available. The best-fit CRM profiles are then assigned to PR pixels to reconstruct the spatial structure of precipitation (yellows in Fig. 1). Figure 5 shows the near-surface rain rate that is estimated in this manner, along with the PR operational product (2A25 surface rain) as a reference. These two retrievals exhibit an excellent qualitative agreement, although some quantitative differences are observed. Comparison of the surface rain rates will be discussed in section 3 in terms of daily statistics.
The nonraining parameters (SST, WV, CLW, and SW) are assigned from TMI brightness temperatures for rain-free footprints that are identified by collocated PR pixels (the algorithm components concerning this part are designated by blue-colored panels in Fig. 1). Figure 6a shows a snapshot of the CWV retrieved for rain-free TMI footprints, where raining TMI footprints are temporarily set as missing. The CRM-provided relative humidity associated with PR-retrieved precipitation profiles defines CWV on raining PR pixels (Fig. 6b). Combining and interpolating those two CWV estimates completes the map of CWV assigned on all PR pixels, as shown in Fig. 6c. The CLW is derived in the same manner, while SST and SW are interpolated only from TMI footprints. It is noted that CWV can be slightly larger in a rain-free area than in an adjacent raining area (e.g., reddish portion in Fig. 6a). This inconsistency can be avoided in the future once an internal retrieval scheme of the nonraining parameters is used, instead of the external data archive that is used in the current study.
The nonraining and raining parameters together are used to simulate microwave brightness temperatures (greens in Fig. 1). Figure 7 depicts observed and computed brightness temperatures (solid and dashed curves, respectively, in the lower four panels), as well as the PR-retrieved raining parameters and assumed DSD and ice-particle models (top four panels) in the nadir cross section for the same scene as Figs. 5 and 6. The ice characteristics will be allowed to be modified in section 3 by varying the ice-particle density ρi, while N0,i is fixed in Eqs. (12) and (13). A displacement of ρi modifies the scattering efficiency of snow and graupel through a corresponding change in both the effective dielectric constant and the slope (λi) in the particle size distribution. A change in ρi is specified by a nondimensional factor fice, or the ratio of ρi to its original value (given in section 2c). The ice-particle model is not yet adjusted in Fig. 7 and, therefore, fice is constant at unity. Figure 7 shows that the computed brightness temperature agrees qualitatively with TMI measurements. A closer examination, however, reveals that the computed brightness temperatures show an underestimation at the low-frequency channels for moderately raining areas (e.g., between scan numbers 1470 and 1500), and an overestimation at the 85-GHz channel over heavily raining areas (e.g., around scan number 1530). These disagreements suggest that the PR-derived precipitation water contents are not sufficiently large to account for either liquid water emissions or ice-water scattering in TMI brightness temperatures in this particular case. The initially assumed DSD model and ice-particle model should be updated to achieve better physical consistency between the PR and TMI. This final portion of the algorithm that adjusts the DSD model and iceparticle characterization is presented in the next section.
3. Adjustment of radar retrieval assumptions
The method introduced in the previous section is extended with an iterative loop to optimize the initial assumption of the DSD model and ice-particle characterization. A simple technique to adjust the initial DSD and ice-particle models is tested here to demonstrate if a reasonable change in these initial assumptions has a sufficient impact on microwave brightness temperatures. Candidate DSD models are the same as those used for the PIA adjustment (Fig. 2). Five different values of ρi are assumed for snow and graupel, where fice (see the last paragraph of section 2g for definition) is allowed to be 1/2, 2/3, 1, 3/2, and 2. The choice of five discrete models of D0 (see section 2c) and fice is rather arbitrary and can be modified in the future to produce a more continuous solution.
The iterative loop starts with the interpolation of microwave brightness temperatures to the PR resolution. This procedure assigns the computed and measured brightness temperatures to every PR pixel. A next guess is assigned to D0 ( fice) in each raining column so that the difference is reduced in Tb,19V (Tb,85V), or the vertically polarized brightness temperature at 19 GHz (85 GHz). The values of D0 and fice at the ith iteration, D(i)0 and f (i)ice, are determined by
where δD0 and δfice are obtained through a Newton–Raphson technique by solving
where T obsb is the observed brightness temperature and T simb is the computed brightness temperature under the assumptions of D(i−1)0 and f (i−1)ice. While ∂Tb,19V/∂D0 is difficult to estimate quantitatively, its sign is always negative because an increase of D0 results in a decrease of the retrieved rainwater content for a given radar reflectivity, which ends up in a smaller emission signal at 19 GHz (and vice versa). The ith guess of D0 is, therefore, determined by decreasing (increasing) D0 to the adjacent value if T simb,19V is colder (warmer) than T obsb,19V. Similarly, a stepwise decrease (increase) of fice defines its ith guess if T simb,85V is colder (warmer) than T obsb,85V, because a lower (higher) ice-particle density produces a weaker (stronger) scattering signal at 85 GHz. The PR matching is performed with an updated assumptions of D(i)0 and f (i)ice, and then brightness temperatures are computed again. This loop is iterated until the brightness temperature difference is minimized at 19 (V) and 85 (V) GHz on all of the raining PR pixels. The minimization of Tb,19V is considered to be achieved when the sign of (T obsb,19V − T simb,19V) is reversed by shifting D0, due to the monotonicity of T simb,19V as a function of D0. The same is true for Tb,85 in terms of fice.
Computational time for the entire process costs 10– 20 times as much as the PR-only algorithm (i.e., only the yellow-colored part in Fig. 1). This ratio could be higher when raining portions within an orbit are small because the “overhead” required for brightness temperature computation over nonraining areas becomes relatively greater in such cases. While it is computationally intensive, the whole process is designed for the purpose of database generation and, therefore, does not affect the computational speed of the TMI rainfall algorithm that will use the generated a priori database.
The above procedure implicitly presumes that a brightness temperature difference at 19 (85) GHz can be attributed primarily to the uncertainty in D0 and fice. This simplification involves a few different problems as follows. First, other microphysical parameters are assumed to be fixed, such as μ in Eq. (1) and N0,i in Eq. (12). The use of an exponential form for ice-particle size distribution has not been justified. Nonsphericity is ignored in the dielectric modeling of hydrometeors. An alternative approach could be to optimize all involved parameters together, assuming that all of the associated uncertainties are given. Such an approach might be better in terms of generality, but would complicate the error propagation through the algorithm, which makes it harder to validate each assumption separately. In practice, the alternative approach would not improve the algorithm performance as long as we have little knowledge to quantify the uncertainties associated with cloud and rain microphysics.
Second, the adjusted assumptions are no longer consistent with the original microphysical assumptions made in CRM simulations. CRM is designed to predict macrophysical variables such as liquid/ice water content (or mixing ratio) with a reasonable accuracy by solving hydrodynamic equations. The treatment of microphysics is much more simplified and less accurate than the macrophysical quantities. It is, therefore, not very meaningful to strictly follow the CRM-provided microphysical characterization, while the macrophysical outputs are more reliable. This fact provides a justification of our approach that D0 and fice are adjusted independently of the original CRM assumptions under a precipitation water content given by CRM as it is.
Finally, uncertainties in the nonraining parameters that are observed within the same TMI FOV as the raining parameters can be also responsible for errors in the computed brightness temperature. A simple diagnosis to justify the assumption adjustment is to see if the brightness temperature difference is reduced also at other frequencies, as well as those used for optimization (i.e., 19 and 85 GHz). If the uncertainty in the DSD is not a major source of errors in brightness temperature, adjusting the DSD would not properly resolve the disagreement simultaneously over different frequencies. Errors in brightness temperatures will be examined later in this section.
Figure 8 shows the retrieved hydrometeor profiles and brightness temperatures for the same satellite snapshot as Fig. 7, but with the adjusted assumptions. The agreement between measured and computed brightness temperatures has improved at every frequency. The computed brightness temperatures at the low-frequency channels have increased because of an increase in precipitation water content for moderately raining regions. This results from the optimization of D0, as seen in the second top panel of Fig. 8, where D0 decreases through the adjustment between scan numbers 1470 and 1490. The initial DSD model is not modified by the TMI adjustment in heavily raining areas (when PIA is reliable). Disagreements between computed and observed brightness temperatures in areas where PIA has been used to adjust DSD would indicate that the discrepancy might be caused by factors unrelated to DSD. Ice scattering signals at the 85-GHz channel around scan number 1530 have deepened in the computed brightness temperature because of an increase in ρi. Although the adjustment of the DSD and ice-density models is performed on the pixel-by-pixel basis, the resultant D0 and ρi exhibit a reasonable spatial variability.
The entire algorithm is also applied to obtain monthly statistics using all of the TRMM orbits during February 1998, with the exclusion of midlatitudinal (i.e., higher than 30°S and 30°N, see section 2d) data. Figure 9 shows the contoured probability density of the correlation between our retrieved and the PR operational (2A25) surface rain rates. Rain rates show excellent correlation when the default DSD and ice-density models (left panel) are used. This simply indicates that the default assumptions, originally taken from the 2A25 algorithm, reproduce the operational 2A25 retrieval. A slight overestimation of our matched CRM solution over the 2A25 rain rate for moderate and heavy rainfall is due to the difference in the treatment of the PIA adjustment between the present and 2A25 algorithms. The contoured probability density spreads out once the assumptions are updated to minimize the discrepancies in brightness temperatures (right panel of Fig. 9). Not only is this expected, it is necessary to account for the systematic underestimation of PR (2A25) relative to TMI in version 5 of the TRMM products.
Figure 10 depicts the computed and observed brightness temperatures (vertical polarization) as a function of the matched CRM rain rate with the updated assumptions of DSD and ice-density characterization. Brightness temperatures increase with rain rate at 10, 19, and 37 GHz until they are saturated, while 85-GHz brightness temperatures decline as rain rate increases. The increase in the low-frequency brightness temperatures is due to the thermal emissions from rainwater, and the decrease at the high-frequency channel is attributed to the scattering by ice particles aloft. The saturation in the emission channels, which occurs at a lower rain rate for a higher frequency, is the combined effect of a large optical thickness of rainfall and the nonuniform beam-filling effect. The observed and computed brightness temperatures show excellent agreement with each other, although the computed 37-GHz brightness temperature has a small positive bias for light rain. This bias arises from a slight inconsistency in the current version of the nonraining parameter retrieval, which results in a small excess in the computed 37-GHz brightness temperatures as mentioned in section 2d. Brightness temperatures with horizontal polarization (not shown) generally exhibit a similar trend to the vertical, except that the discrepancy is consistently larger by a few kelvins for horizontal polarization than for the vertical at 85 GHz. This can be attributed to microwave scattering by nonspherical ice particles, which is not taken into account in the present analysis.
Figure 11 shows rms errors in computed brightness temperatures relative to the measurements. Discrepancies have been reduced by adjusting the DSD and ice-density assumptions at all the TMI channels, despite the fact that D0 and fice were adjusted using only the vertical polarization at 19 and 85 GHz, respectively. As mentioned earlier in this section, the fact that the optimized D0 and fice result in a coherent impact over different frequencies and polarizations is encouraging for the legitimacy of the assumption adjustment approach. Figure 11, however, also exhibits consistent errors remaining after the optimization of D0 and fice, particularly at 85 GHz. The remaining errors are ascribed to other contributions, such as the nonraining parameters and cloud/precipitation microphysical variables, which are not considered here. Refining the retrieval algorithm of the nonraining parameters and the treatment of microphysical parameters is being planned for future improvement.
4. Concluding remarks
This paper introduces a methodology to analyze precipitation profiles based on TRMM PR and TMI measurements, focusing on describing the architecture of the algorithm in detail. The raining parameters are derived from PR reflectivities with assumed DSD and ice hydrometeor models, while the nonraining parameters are retrieved from TMI brightness temperatures for rain-free areas. The PR PIA provides a constraint on the DSD assumption for heavy rainfall. The raining and nonraining parameters together are then employed for computing brightness temperatures to be compared with TMI measurements. Differences between computed and measured brightness temperatures are reduced by iterating the entire procedure with updated assumptions of the DSD and ice-density models.
The current analysis provides a PR-generated database without any TMI adjustments and a combined PR–TMI database following the brightness temperature adjustment, separately from each other. A ground validation effort may be necessary to verify that the TMI-adjusted DSD is indeed closer to observations than the initial guess. The adjusted database, in parallel to the CRM-only database, can serve as the a priori database in the GPROF algorithm (Kummerow et al. 2001). It is expected to represent a significant improvement over the current CRM-only formulation and is likely to go a long way in understanding discrepancies between radar-only and radiometer-only solutions. This will be investigated in a follow-up paper.
This research is supported by NASA’s Precipitation Program Grant NAG5-13694. The TRMM datasets (1B11, 1C21, 2A12, and 2A25) were provided by the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center (GSFC). The nonraining parameter dataset retrieved from TMI was obtained from the Remote Sensing Systems (RSS) Web site (online at http://www.ssmi.com/tmi/tmi_browse.html), sponsored by NASA’s Earth Science Information Partnerships (ESIP)—a federation of information sites for Earth Science—and by NASA’s TRMM Science Team.
Derivation of PR Radar Equation
Radar reflectivity Z is related with radar-received powers Pr through the radar equation for scattering particles,
where Pt is the transmitted power, Gt the transmit antenna gain, Gr the receive antenna gain, h the pulse length, θd the down-track beamwidth, θc the cross-track beamwidth, λ the wavelength, r the distance, and
the scattering property of particles with their dielectric constant ε.
Equation (A1) is designed so that Z is equivalent to the sixth moment of DSD if 1) scattering particles are so small that the Rayleigh approximation applies, 2) ε is independent of the water temperature, and 3) the attenuation is negligible. Raindrops at the PR frequency, however, satisfy none of these conditions, in general (L’Ecuyer and Stephens 2002). Hence, a more accurate form of the radar equation should be introduced to describe the actual relation of DSD with the PR received power,
where σsb,i, ksext,i, b, and kext are the single-particle backscattering and extinction coefficients of the ith hydrometeor species and the total backscattering and extinction coefficients, respectively. The total number of the hydrometeor species nspec is as much as 6 in the present study, which consists of rain, snow, cloud water, cloud ice, graupel, and the melting particles. The Mie approximation is assumed to calculate σsb,i and ksext,i for all the species. The combination of Eqs. (A1) and (A2) leads to
or, for a nadir satellite view,
is the extinction optical thickness between z1 and z2, and
is the correction factor accounting for the self-extinction of a scattered echo within a radar pulse with the length of Δz. Assuming that b and kext are constants within the vertical interval of Δz, one can perform integrals in Eq. (A8) to obtain
The wavelength λ equals 21.7 mm for the PR frequency of 13.8 GHz and |K| of 0.9255, the value for liquid water of 0°C at 13.8 GHz, is assumed in the TRMM 1C21 dataset to evaluate the radar reflectivity. The PIA is calculated as
Corresponding author address: Dr. Hirohiko Masunaga, Department of Atmospheric Science, Colorado State University, Fort Collins, CO 80523. Email: firstname.lastname@example.org