Microphysical Properties of Frozen Particles Inferred from Global Precipitation Measurement ( GPM ) Microwave Imager ( GMI ) Polarimetric Measurements

Scattering differences induced by frozen particle microphysical properties are investigated, using the vertically (V) and horizontally (H) polarized radiances from the Global Precipitation Measurement (GPM) Microwave Imager (GMI) 89 and 10 166 GHz channels. It is the first study on frozen particle microphysical properties on a global scale that uses the dual-frequency microwave polarimetric signals. From the ice cloud scenes identified by the 183.3±3 GHz channel brightness temperature (TB), we find that the scattering by frozen particles is highly polarized with V-H polarimetric differences (PD) being positive throughout the tropics and the winter hemisphere mid-latitude jet regions, including PDs from the GMI 89 and 166 GHz TBs, as well as the PD at 640 GHz from the 15 ER-2 Compact Scanning Submillimeter-wave Imaging Radiometer (CoSSIR) during the TC4 campaign. Large polarization dominantly occurs mostly near convective outflow region (i.e., anvils or stratiform precipitation), while the polarization signal is small inside deep convective cores as well as at the remote cirrus region. Neglecting the polarimetric signal would easily result in as large as 30% error in ice water path retrievals. There is a universal “bell-curve” in the PD TB relationship, where the PD amplitude peaks at ~ 10 K for all three channels in the tropics and increases slightly with latitude (2-4 K). Moreover, the 166 20 GHz PD tends to increase in the case where a melting layer is beneath the frozen particles aloft in the atmosphere, while 89 GHz PD is less sensitive than 166 GHz to the melting layer. This property creates a unique PD feature for the identification of the melting layer and stratiform rain with passive sensors. Horizontally oriented non-spherical frozen particles are thought to produce the observed PD because of different ice scattering properties in the V and H polarizations. On the other hand, turbulent mixing within deep convective cores inevitably 25 promotes the random orientation of these particles, a mechanism works effectively on reducing the PD. The current GMI polarimetric measurements themselves cannot fully disentangle the possible mechanisms.

and 166 GHz (right panels) polarimetric channels. The GMI observed radiance distributions and polarimetric differences (PDs) are shown in the top and bottom, respectively.
Capable of penetrating clouds though, the PDs from 89 and 166 GHz contain surface polarization signals that need to be 5 taken into account prior to cloud analyses. The polarized surface contributions are particularly strong over water where the degree of linear polarization (DoLP) of microwave radiances depends on surface roughness induced by wind and sea foams (Meissner andWentz, 2012, Chen, 2014). As seen in Fig. 1c, the 89 GHz PDs exhibit a large contrast between land and ocean, but such a contrast disappears at 166 GHz. This feature is again manifested in Fig. 2, where observations over ocean, marked by black open circles, spread to large PDs at the warmest end from 89 GHz observations, while this branch is not distinguishable at 10 166 GHz.
To study scattering properties of cloud and frozen particles jointly with the 89 and 166 GHz, we need to isolate cloud cases from the clear sky where the surface contributions are not negligible. To flag heavy cloudy cases, we use the simultaneous 183.3±3 GHz measurements, which show little sensitivity to surface emission due to the strong water vapor absorption but still with enough sensitivity to deep convective clouds. As described in Gong and Wu (2013), a "3σ method" can be used to distinguish 15 between clear-and cloudy-sky scenes for the Microwave Humidity Sounder (MHS) tropical observations. The "3σ method" first identifies the peak (TB 0 ) and standard deviation (σ) of the Probability Density Function (PDF) of a month-long radiance Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-787, 2016 Manuscript under review for journal Atmos. Chem. Phys. Published: 21 September 2016 c Author(s) 2016. CC-BY 3.0 License. observation in a given region or a latitude belt, and then the cloudy-sky scene is defined at places where TB < TB 0 -3σ. In this study, we apply the same methodology to select cloudy-sky scenes. Cloud occurring frequency after applying this "3σ method" is comparable to that observed by CloudSat radar (see Appendix D for a detailed map). Because of ocean/land differences, separate PDFs are computed for ocean and land surfaces in each 10° latitude bin. The three GMI channels (89, 166 and 183 GHz) have roughly the same (7.2 km × 4.4 km) footprint size, which minimizes the sampling error due to beam-filling effects 5 between different frequencies. The combined sensitivities of the three channels will help to distinguish the scattering differences between frozen precipitation (large particles) and ice clouds (small ice crystals).
As a caveat, GMI swaths are mismatched between the LF and HF channels, offset by as much as 4 degrees in pointing (L1C, ATBD). Therefore, the L1C data are not appropriate for joint analyses of 89, 166 and 183 GHz measurements. Instead, we use the L1C-R product that is derived from the L1B radiances using the nearest-neighborhood matching method to register 10 collocated LF and HF observations. This re-registration procedure eliminates ~ 7% of the data, of which most occur at the edge of the LF scan swath. These missing data do not impact our climatological study. A comparison of the L1C-R and L1B datasets is given in the Appendix A.
While the 89 and 166 GHz polarimetric measurements are useful to infer microphysical properties of frozen precipitation, PD at a higher frequency can reveal cloud properties with smaller ice crystals. To expand our study to ice clouds with small ice 15 crystals (e.g., cirrus), we also include an analysis of the 640 GHz V-and H-pol data acquired by NASA airborne Compact Scanning Submillimeter-wave Imaging Radiometer (CoSSIR) during the Tropical Composition, Cloud and Climate Coupling (TC4) campaign in July-August 2007 (Evans et al., 2005;. The CoSSIR instrument had a similar conical scan geometry to GMI in both forward and backward directions. About ~39000 CoSSIR conical samples from three TC4 flights are analyzed.
These samples were mostly over the Pacific Ocean near Central America and large enough to commensurate with the GMI 20 statistics over the tropical oceans. We found that CoSSIR measurements were unreliable during the ascent/descent flight in altitude, and therefore only measurements at a cruise altitude above 19 km are used.

Data analysis
PD-TB v relationship is a very useful diagnostic tool and will be used throughout this study. Fig. 2 shows a scatter plot of PD vs. TB v from the same squall line event reported in Fig. 1. As expected, for large oceanic PDs seen at 89 GHz, there is a 25 distinct branch in Fig. 2a, showing a strong linear PD-TB v relationship at the warm or clear-sky TBs. This branch is separated from the frozen particle scattering cluster that exhibits a bell curve in the PD-TB v relationship. The bell curve is characterized by a peak of PD = ~10 K at TB v = 220 K. Note that the surface branch is absent at 166 GHz due to strong water vapor absorption in the subtropics (Fig. 2b), but the bell curve is similar to the one at 89 GHz with peak of PD = ~10 K at TB v = 190 K. The common bell curve indicates that PD is small at the warmest (i.e., thin cloud or clear-sky) and coldest TB (i.e., deep convection) cases. As 30 shown in the later sections, the bell curve appears to be universal at all latitudes.
Atmos. Chem. Phys. Discuss., doi:10.5194/acp-2016-787, 2016 Manuscript under review for journal Atmos. The PD climatology is computed globally and sorted according to latitude, TB v , and PD value itself in 10°, 10K, and 1K bins, respectively. We also compiled a two-dimensional (2D) PDF with respect to PD and TB v for each latitude bin. Monthly statistics of the PD -TB v relationship are characterized by the peak and standard deviation of each 2D PDF. Only the cloudy-sky portion of the PD -TB v relationship is reported in this study using the 3σ cloud detection with the 183.3±3 GHz radiances in each 20 latitude band and computed separately for land and ocean.
Six months of GMI L1C-R data (June, July, and August in 2014 and 2015) have been compiled to study the cloudy-sky PD statistics. Although we use an earlier data version (V03) in this study, we find that the results are similar between V03 and V04 and choose to use the V03 data for the rest of paper. For the CoSSIR 640-GHz PD data, of which the campaign was based in Costa Rica, the statistics will be compared to the GMI data from the [0, 10°N] latitudinal bin.

25
To study the scattering properties of cloud ice and precipitation-sized frozen particles at 89 and 166 GHz, we make use of precipitation measurements from GPM DPR in our analysis. Galligani et al. [2013] found that large 85 GHz PDs of TRMM-TMI were often associated with a melting layer underneath (detected by TRMM-PR), suggesting a strong connection between cloud ice processes and surface precipitation. Thus, we employ the GPM Ku-band radar bright band (BB) flag to evaluate the 89 and 166 GHz PD properties and their connection to the presence of BB, as BB is always associated with stratiform precipitation. For this analysis, we first collocate the measurements between the Ku-band radar and the GMI 89 and 166 GHz data, and then tag the BB flag to each GMI measurement. The Ka-band radar also produces a product for classifying stratiform and non-stratiform precipitation cases, but the Ku radar BB flag is generally considered to be more reliable (Iguchi et al., 2010).

5
To evaluate impacts of frozen particle orientation and ice water content (IWC) on the TB and PD -TB v relationship, two numerical RTMs are employed in this study. The first model, called RT4, is a variant of the Polarized Radiative Transfer Model Distribution (PolRadTran) developed by Evans and Stephens (1995). For our study, this model is configured to mimic GMI and CoSSIR's viewing geometry and channel frequencies, as well as to include several ice habitats (snowflake dendrites, hexagonal plate, cylindrical columns, bullet-rosettes with 7 branches) and their scattering properties for diameter ranging from 20 to 2000 10 µm as described in Yang et al. [2013]. As described in Evans and Stephens (1995), the discrete dipole approximation (DDA) method is used to speed up the scattering calculation once the ice refractive index is provided. This index has weak dependence on temperature and we used the value at 240K in this study. Values at other temperatures (220K, 260K) have also been tested, and only show minor differences in the results. Frozen particle obey a Gamma size distribution, which is characterized by mean mass-weighted equivalent sphere diameter (D me ). The desired gamma size distribution is achieved using an optimization 15 procedure to adjust the equivalent volume diameter (D e ), where "# = # % / # ' . This optimization procedure assures that the truncated and discrete size distribution has the desired moments (i.e. D me and dispersion width).
The RT4  The RT4 model has several limitations. The modeled ice crystal shapes (except dendrites) represent primarily the pristine mid-latitude cirrus from synoptic outflows, which is not ideal for larger ice particles in the stratiform precipitating regions nor for tropical deep convection cases. While the DDA calculation is very efficient for computing cirrus particle scattering, it is still 25 computationally expensive for large particles, of which the valid range of IWP is limited to <1 kg/m 2 and D me to <500 µm.
Besides, the RT4 surface model and homogeneous ice cloud layers are often over-simplified compared to the reality. Because of these limitations, RT4 is perhaps mostly suitable for simulating the 640 GHz PDs where the surface contributes the least, compared to the lower frequency channels, and this channel is insensitive to precipitation-sized ice particles. Assuming randomly particle orientation, Evans et al. [2012] applied the RT4 model with a dynamical mixture of different particle shapes to cloud ice 30 retrievals on the CoSSIR data.
Another RTM for simulating frozen particle scattering, called Cloud Radiance Model (CRM), is based on the cloudy-sky forward model in Aura Microwave Limb Sounder (MLS) (Wu and Jiang, 2004;Gong and Wu, 2014). CRM is a non-polarized model, and it assumes only spherical ice particles in the scattering calculation. It has been applied to simulate cloud-induced radiances at frequencies < 3 THz for a complicated viewing geometry (e.g., satellite nadir/limb sounding, in-cloud aircraft flight, 35 and ground stations). Like RT4, CRM assumes the Gamma size distribution in the current study. Different from RT4, CRM allows the user to specify the vertical distribution of IWC or liquid water content (LWC) in the mixed-phase case, to study sensitivity of the cloud-induced radiance to cloud layer height. Cloud height is another free parameter in addition to particle size and IWP, which can alter the cloud radiance amplitude significantly in similar frequency range   10 This is also evident at 89 and 166 GHz over land (Figs. 4c,d), which may be associated with the drier atmosphere at highlatitude in the winter atmosphere. Although the dynamic range of TB V shrinks with increasing latitude, the PD peaks are still detectable and occur at about the same TB V value (~220K) for the 89 GHz oceanic clouds credited to the large clear-sky variability and the resulting 3σ threshold. These findings suggest that the ensemble characteristics of the ice crystal shape and orientation that 89 GHz is sensitive to are consistent across all latitudes, which is a piece of good news to both the satellite 15 retrieval and cloud modeling group, as many cloud microphysics assumptions during the retrieval and modeling processes are based on mid-latitude observations. The stringency of shape and orientation parameters ensures less uncertainty in the retrieving/modeling of the tropical and high-latitude atmosphere.  climatologies. Top panels for ocean and bottom panels for land conditions, respectively. Surface signal is removed using the 183.3±3 GHz radiance and the "3s threshold" method.
The oceanic 166 GHz PDs at TB V = 200-220 K appear to increase with latitude in both hemispheres (Fig. 4b). There is also 5 a hint in Fig. 4a for 89 GHz although the tendency is rather weak. This latitudinal dependence of the PD amplitude could be due to latitude-dependent water vapor abundance, because ample water vapor can damp the PD signal. The water vapor damping is stronger in the tropics and consequently flattens the PD-TB curve more. The water vapor attenuation effect is also evident for ocean conditions at mid-high latitude of the winter atmosphere (30°S -70°S), which is caused by the fact that the atmosphere there is very dry and cannot block the surface signal anymore (e.g., Fig. C1d). However, the upper troposphere water vapor 10 amount should be roughly symmetric between the two hemispheres (Solden and Lanzante, 1996), but the PD curves do not suggest such an inter-hemisphere symmetry. Over ocean, the PD peak amplitude is slightly larger in the summer hemisphere, while is contrary over land. Therefore, the role of water vapor in attenuating the PD is perhaps more complicated than expected.
We will explore this in more detail in the next section.  15 of frozen particles. In the next section, we will discuss two plausible mechanisms that may explain the observed features.

A Simple Model for the PD-TB relationship
The leading hypothesis of the observed positive PD in cloudy sky is that the bulk ice properties produce different volume  remote sensing techniques. Clouds may produce little or no PD if the shaped particles are randomly oriented (e.g., tumbling in a turbulent environment). On the other hand, as discussed below, the PD tends to decrease with cloud optical thickness as a result of cloud self-extinction in the radiative transfer.
Let us assume that frozen particles are horizontally aligned, and the V-and H-pol extinction coefficients are dominated by scattering. Naturally V-and H-pol radiances pass through the same cloud with slightly different optical depths (τ V and τ H ). We 5 define the aspect ratio (AR) factor as ≡ , -to characterize this optical depth difference. Fig. 6 illustrates a simple twolayer model, where cloud optical depth is τ 2 . Neglecting the water vapor attenuation, radiance at the top-of-atmosphere (TOA) can be expressed analytically as,

25
This simple 2-layer model (Eqn. 3) is able to reproduce several fundamental features of the observed PD-TB relationship.
As shown in Fig.7a Because T J -T 1 is roughly constant, this linear dependence on τ 2V -τ 2H helps to spread the slopes for different ARs because τ 2V - The model can also explain the complex PD-TB relationships seen in the GMI oceanic data where the background surface emission is polarized. Because the 89-GHz channel is more sensitive to the polarized surface emission than the 166-GHz, the observed cloud PD often contains both surface and frozen particle contributions. As shown in Fig. 7b, the combined PD from the surface and cloud scattering can result in a tilted bell curve. The observed PD statistics in Fig. 3 are likely a mixture of the two cases shown in Fig. 7.

10
Finally, the conceptual model predicts weak dependence on channel frequency for the PD-TB relationship. The reason is as follows. If cloud extinction is dominated by scattering, τ is proportional to the 4 th order of the frequency, i.e., = % . Therefore AR is independent of frequency by definition. If the background emission is unpolarized (i.e., 1-C = 1, C = 1 ), the PD from Eqns. 1-2 has a maximum of ( 1 − @ )( 34 56 − 3DE * 4 56 ) at @-= ln ( ) ( − 1), which is only weakly dependent on frequency. This weak-dependence of PD on frequency seems to be consistent with the observations that PD amplitudes barely 15 increase from 89, 166 to 640 GHz.  Although AR in the range of [1.2, 1.4] produce a good envelope for the observed bell curves, AR = 1.3 seems to fit best overall 5 to the observations. This AR value was also found in Davis et al. (2005) in studying the 122 GHz data. Having said that, there is an indication that the best fit of AR might increase with frequency, suggesting an increasing contribution of horizontally-oriented particles at higher frequencies. Since the higher frequencies are more sensitive to small frozen particles, it implies that small frozen particles have an increased percentage of the horizontal orientation. While it is difficult to infer the detailed composition of particle habits in these clouds, ice crystals in cirrus clouds (more sensitive by 640 GHz) are likely to have a predominant type 10 (e.g., plate, column, etc.) and snow aggregates (more sensitive by 89 GHz) tend to have diverse irregular shapes. It is also found in the simulation that the modeled PD maximum decreases slightly at 640 GHz compared its value at 166 GHz, which is indicative of stronger water vapor attenuation at 640 GHz.

Full Polarized Model Simulation
AR is a parameter that characterizes the collective scattering property of frozen particle habitats. To study the relative 15 importance of individual particle habitat, we use a fully polarized RT4, a variant of the PolTranRad model [Evans and Stephens, 1995], to simulate the PD-TB relationship. Four types of ice habitats and different D me are used in the simulations, and the results are shown in Fig. 9 for 640 GHz assuming all frozen particles are horizontally aligned. As seen in Fig. 9, snowflake dendrites and cylindrical columns are most effective shapes to produce large PD values at 640 GHz, whereas cylindrical columns and bullet-rosettes can also induce a similar PD amplitude to what is measured by CoSSIR. The first two ice crystal types result in 20 excessively large PDs in the PD-TB relationship. For bullet-rosettes, on the other hand, the PD does not produce the bell curve when the cloud becomes optically thick, whereas the cylindrical column shaped crystal with D me between 200 and 300 µm produces the best agreement with the observation. According to Yang et al. [2013], this type of ice crystal with D me = 250 µm, when projected to GMI's viewing direction, has an AR of 1.46, a value in a surprising good agreement with that from the CRM simulations shown in Fig. 8c. This consistency also supports the AR values used in the simple conceptual model and its radiative 25 transfer argument for the bell-shape curve. We hope with more accurate multi-frequency polarimetric measurements that ice habitats may be retrievable from the PD -TB relationship eventually. More questions also rise up. For example, given that ice  In the rest panels, red, yellow, cyan, and blue denote D me = 200, 300, 400 and 500 µm, respectively.
One of the top questions investigated in this study is: what process causes the low PDs of the bell curve at cold TBs? One would argue that at cold TBs (i.e., strongly depressed by the scattering of deep convective clouds) the orientation of frozen 10 particles were mostly random instead of horizontally aligned. As a result, those deep convective cold TB cases would yield a weaker PD. PolTranRad has another version called RT3, which allows to simulate effects from randomly orientated ice crystals, and the bell curve is found to be still present, despite the amplitude is rather small (< 3K) regardless of type of ice habitat used.
Thus, the random orientation of frozen particles does not completely cancel out the particle shape effects, as revealed by the RT3 simulations.

15
A varying mixture of randomly and horizontally-oriented particles may also provide a plausible explanation of the observed bell curve. If the larger mixing factor "α" corresponds to a more mixed environment inside clouds, α is likely to increase with deep convection strength, or a colder TB. On the other hand, horizontally aligned non-spherical particles are likely more effective to produce PDs in the anvil outflows than convective cores, because random orientation within the turbulent convective core tends to cancel out the PD. This mechanism seems to have a support evidence in the ground-based Doppler radar 20 observations (Homeyer and Kumjian, 2015), where they often found lower differential reflectivity (ZDR) near the deep convective core but larger values in the anvil outflows. In addition, the irregular shape of graupel found near convective cores likely contributes to a small PD as well as suggested in their paper.

PDs from melting layers and ice clouds
Since 89 and 166 GHz radiances are sensitive to scattering and emission of the precipitation-sized frozen particles, we need 25 also to take into account the effects of the floating snow layer (Xie et al., 2012) or the melting layer close to the surface (Galligani et al., 2013) underneath the cloud layer. Although it is difficult to quantify the PD sensitivity separately for different layers of frozen particles, it is possible to evaluate the response of PD to the melting layer with the combined DPR-GMI observations. Melting layers manifest themselves as a unique bright band (BB) on the Ku-band radar reflectivity profile, as a result of stratiform rainfall. We therefore use the BB flag and precipitation flag in DPR's Ku-band radar Level-2 products to discriminate the stratiform and non-stratiform precipitation scenes. Non-stratiform precipitation types include convective rainfall and snowfall.  should be related to lightly precipitating scenes that otherwise would show much stronger PD signals. In the case with BB, most of the polarized surface emission is blocked by BB. For non-stratiform precipitation (without BB), most likely there will be no PD signal from both channels (i.e., the center of the black contours is near the origin point in Fig. 10a). This indicates that the dominant PD signals of 89 and 166 GHz indeed come from ice cloud for non-stratiform precipitation, even it is snowfall. When 15 the melting layer is present, both 89 and 166 GHz PDs incline to be positive. As a result, when the melting layer is present, the peak of the joint PDF (dark red) moves to ~ 5K and 8K for 89 and 166 GHz, respectively, as shown by the color shades in Fig.   10a, implying that the melting layer in fact contributes partially to the PD signal. This finding is consistent with Galligani et al. (2013) where they found that 85 GHz PD of TRMM's TMI is likely positive when the melting layer is present. It is also consistent with the higher BB ZDR as observed by ground-based radars (e.g., Oue et al., 2014). While 166 GHz is more sensitive 20 to smaller ice crystals than 89 GHz, the more positive PD response in the 166 GHz indicates that smaller ice crystals should be more obsolete and/or more horizontally oriented. Larger particles like snow aggregates likely exhibit irregular shapes and the scattering signal in 89 GHz would hence be less directional, resulting a smaller PD response in this channel. This finding is highly consistent with the weak increase of AR value with the increase of channel frequency in CRM simulations shown in Section 4.1. Moreover, frozen particles of all sizes tend to be more horizontally oriented for stratiform rainfall compared with Over land, the PD characteristics vary roughly the same way with that over ocean, but the 4:1 branch exhibited in the ocean case is not present in the land non-stratiform precipitation cases. The PDF peaks along the 1:2 line are different in the cases with and without BB. Other factors, e.g., liquid emission, may also play a non-trivial role here, but they are hard to be justified at this moment. Due to the complexity of melting layer, this paper does not intend to further quantify its role on the PD's characteristics.

Conclusion and Future Work
High-frequency microwave polarimetric measurements contain ample information of frozen particle habitat and orientation.
In this paper, the GMI 89 and 166 GHz PD characteristics are analyzed, which reveal a universal bell-curve in the PD -TB relationship. The bell curve also exists in the CoSSIR 640 GHz observations. All three frequency channels show a maximum PD amplitude of ~ 10 K, occurring around TB = 200 -220 K. Such a PD-TB relationship is robust for almost all frequencies and 10 latitudes. Moreover, both 89 and 166 GHz channels show positive PD responses when a melting layer is present. This is a new feature not reported before, suggesting that at least part of the PD signals comes from the melting layer. The melting layer contribution is confirmed by comparing the PD statistics with non-stratiform precipitation situations.
Using a simple analytical radiative transfer model, we can explain the bell curve of the PD -TB relationship with the term of aspect ratio (AR). In this simple model we assume that the V-and H-pol extinction coefficient profiles are similar and scaled by dataset that does not collocate the LF and HF channels, 183.3±3 GHz TB now corresponds much better with 89 GHz TB.
Moreover, 166 and 89 GHz TB shows a much better coherence in the L1C-R data, showing two separate branches more clearly in Fig. A1 (d). The different slopes of these two branches correspond to the clear-sky and cloudy-sky characteristics, 10 respectively, and the clear-sky branch also is the warmest in 183.3±3 GHz TB. Scatter behavior of Fig. A1(d) is very similar to comparable instrument, e.g., Special Sensor Microwave Imager and Sounder's 91 and 150 GHz channels (not shown).

Appendix B: Estimation of Instrument Noise and Distribution of Negative PD
GMI's 89 and 166 GHz PDs are not always positive. As a matter of fact, negative values occur frequently, and their occurrence is independent with TB. However, large negative PD only happens occasionally. The histogram of negative PD 15 closely follows the Gaussian distribution, as shown in Fig. B1 by the black solid fitted lines. The "Gaussian_fit" function from IDL software package is used for this purpose, where the mean (µ) and standard deviation (s) parameters are given as output.
The standard deviation of the fitted Gaussian curve is a good estimate of the instrument noise (e.g., Wu et al. [2009] If the negative PDs are purely originated from instrument noise, they should be randomly distributed geographically. However, that is not the case for 166 GHz. If we plot out the occurrence frequency (OF) map of negative PDs that are smaller than the negative instrument noise estimated using the aforementioned method, the OF is trivial at 89 GHz, indicating that 5 negative PD is likely instrument noise. But for 166 GHz, large negative PDs occur quite frequently (up to 18% of the time), and they seem to favor the winter-hemisphere side of the tropical deep convections, as shown by Fig. B2. This is the first time such a geographical preference has ever been reported.
One theory explains the negative PD as the scattering effect from predominantly vertically oriented ice crystals that may be oriented by the electromagnetic field created by lightning near the cloud top (Prigent et al., 2005;Homeyer and Kumjian, 2015).

10
We can also consider the vertical orientation a counterpart solution of the horizontal orientation from the RTMs as the latter can only produce positive PD signals. However, global lightning distribution, which essentially connects to the tropical deep convective activities, correlates poorly with the geographic distribution of negative PDs (e.g., Cecil et al., 2014). Instrument noise is not likely the underlying cause for 166 GHz either since it would be randomly distributed otherwise. Having said that, the secondary branch in the negative PD branch in Fig. 4 surely raised up more interests on studying the causes of negative PD.

15
Although the actual cause remains unclear at this moment, negative PD, especially the large magnitude one, assures a rather attractive topic for future studies.  method is proven here to be able to effectively exclude the surface signal in the PD -TB relationship, although it might be too stringent, especially for 166 GHz, as the clear-sky PD is only weakly dependent on the ocean surface roughness for this channel.

Appendix D: Cloud Occurring Frequency (CF) derived from the "3σ method"
To demonstrate the credibility of the "3σ method" on screening out the clear-sky scenes, we show the CF of July 2015 from the 89 and 166 GHz channels in Fig. D1. The geographic distribution of CF of anvil and deep convective clouds agrees well with 5 that observed by the CloudSat radar (Sassen et al., 2009). Due to the stringent criteria set by the "3σ method", many cirrus scenes are classified as clear-sky scenes, and therefore cirrus cloud CF is obviously biased low using this method.