Quantitative Investigation of Radiometric Interactions between Snowfall, Snow Cover, and Cloud Liquid Water over Land

Falling snow alters its own microwave signatures when it begins to accumulate on the ground, making retrieval of snowfall challenging. This paper investigates the effects of snow-cover depth and cloud liquid water content on microwave signatures of terrestrial snowfall using reanalysis data and multi-annual observations by the Global Precipitation Measurement (GPM) core satellite with particular emphasis on the 89 and 166 GHz channels. It is found that over shallow snow cover (snow water equivalent (SWE) ≤ 100 kg m−2) and low values of cloud liquid water path (LWP 100–150 g m−2), the scattering of light snowfall (intensities ≤ 0.5 mm h−1) is detectable only at frequency 166 GHz, while for higher snowfall rates, the signal can also be detected at 89 GHz. However, when SWE exceeds 200 kg m−2 and the LWP is greater than 100–150 g m−2, the emission from the increased liquid water content in snowing clouds becomes the only surrogate microwave signal of snowfall that is stronger at frequency 89 than 166 GHz. The results also reveal that over high latitudes above 60°N where the SWE is greater than 200 kg m−2 and LWP is lower than 100–150 g m−2, the snowfall microwave signal could not be detected with GPM without considering a priori data about SWE and LWP. Our findings provide quantitative insights for improving retrieval of snowfall in particular over snow-covered terrain.


Introduction
Passive microwave (PMW) retrieval of snowfall is one of the most challenging components of precipitation monitoring from space, with the largest error in precipitation retrieval often related to snowfall [1][2][3][4][5][6] over snow cover [7]. Snowfall emission is almost negligible due to the low dielectric constant of ice particles, especially over emissive land surfaces. Falling snow and ice particles scatter the upwelling surface radiation at high microwave frequencies and thus decrease the observed brightness temperatures (Tb) at the top of the atmosphere. This radiometric signal, however, is much weaker than the overland rainfall scattering [8][9][10] and can be significantly masked due to the confounding effects of increased cloud liquid water path (LWP) during snowfall and reduced surface emissivity as a result of snow accumulation on the ground.
Snow particles have complex and variable nonspherical shapes and bulk densities. The size distribution of snowfall particles depends on numerous factors, including temperature, pressure, and the level of water vapor path (WVP) supersaturation at different vertical layers of the atmosphere [11][12][13][14][15]. The nonspherical snow particles usually have lower densities than raindrops with equal mass, which causes them to exhibit weaker scattering [9,10,16].
Radiative transfer modeling shows that this weak scattering becomes detectable at frequencies above 80 GHz and reaches its maximum at frequencies 150 to 166 GHz [17,18]. However, in snowy clouds thermal emission from the presence of supercooled liquid water content can mask the already weak snowfall signal even at these high-frequency channels [19]. The warming effects of the liquid water emission in the mixed-phase clouds may even exceed the cooling effects of snowfall scattering [20] and completely mask the snowfall signal [21,22], which adds to the complexity of (light) snowfall retrievals.
The presence of snow cover on the ground is another challenge in detecting the snowfall scattering signal. Snow cover is a relatively strong scatterer at microwave frequencies above 20 GHz [23,24], and this scattering increases monotonically with frequency up to 100 GHz [25]. Increased scattering reduces the surface emissivity depending on physical and microphysical properties of snowpack such as depth, density, wetness and the distribution of grain size [24,26,27], which vary in response to the snow metamorphism [28][29][30][31]. A thicker and denser snowpack often scatters the upwelling surface emission more, especially when snow ages and develops larger and denser particles [23,27,32,33]. In addition, radiometric properties of snow cover are very sensitive to its liquid water content. For a very small amount of liquid water content of around 2%, the absorption dominates scattering and turns the snowpack to almost a blackbody radiator [34]. The combination of the explained radiometric processes has two important consequences, which add to the complexity of PMW snowfall retrievals. First, there is a likelihood that the snow cover and snowfall microwave signatures become very similar [2,35,36]. Second, the low emissivity of dry snow cover can significantly weaken the already weak snowfall scattering [6,[37][38][39] or vanish it completely.
Radiative transfer modeling of both snowfall and snow-cover has large complexity under diverse environmental conditions over land with large variability over temporal and spatial scales and requires many input parameters [40]. Most radiative transfer modeling has been conducted assuming that the snowfall particles are spherical [1,15], and single scattering theories can approximate their scattering. The discrete-dipole approximation is also used to account for the nonspherical shape of snowfall ice particles [41,42]. Nevertheless, snowfall radiative transfer models can often account only for a limited number of snowfall particle shapes [43,44] and lack the ability to properly address bulk scattering of a snowfall profile throughout the atmospheric column [45]. Additionally, modeling of the emission signal of the supercooled liquid water content in snowing clouds is still not well parameterized, especially at frequencies above 31 GHz [15,46]. Using the Rayleigh approximation in the absence of precipitation for water droplets larger than 0.2 mm, the authors of [47] calculated the microwave absorption of supercooled liquid water at 21 and 31 GHz and found that this absorption strongly depends on cloud temperatures. Thus, the resulting absorption derived from common dielectric models significantly deviates during the snowfall at temperatures below 270 K mainly due to poor representations of the primary relaxation frequency of water.
The snow cover scattering below 100 GHz can be explained by the Born approximation approach [48], which partially accounts for the near-field radiation produced by adjacent snow grains, using the low-frequency practical media theory [49]. More recently, in order to extend the range to larger particles and higher-frequencies, Grody [27] proposed the use of a quasi-crystalline approximation to bridge the gap between small-and large-size snow grains. However, the radiative transfer modeling of snow cover effects is not conducted for frequencies above 150 GHz [30] in which the primary snowfall retrieval bands are located in the PMW spectra. The uncertainty of this high-frequency approximation is often quite larger than that at lower frequencies because the optical parameters of snow cover become less frequency-dependent and saturated as the particle diameter approaches the wavelength.
Clearly, complexities of radiative transfer modeling are amplified when it comes to modeling snowfall over snow cover, accounting for the presence of cloud liquid water content, where there is significant heterogeneity in the characterization of the initial and boundary conditions at a global scale. The Global Precipitation Measurement (GPM) Mission PMW observations can help to unpack these relatively unknown radiometric interactions and to provide some important insights into PMW snowfall retrievals. A recent empirical study conducted by [38] using coincidental CloudSat and GPM observations found that cloud liquid water emission increases the brightness temperatures up to 10 K at 166 GHz and usually even more at 89 GHz. In addition to these interesting findings, previous studies have not considered the synergistic effect of atmospheric LWP and surfaceaccumulated Snow Water Equivalent (SWE) on the snowfall signal.
The goal of this paper is to quantify the radiometric effects of snow cover and cloud liquid water content on the microwave signatures of snowfall at high-frequency channels, largely focusing on observations from the GPM core satellite. The findings of this paper push the understanding of the land-atmospheric effects one step further, investigating the conditions in which the passive radiometric signal of snowfall is affected by surface and atmospheric components and thus cannot be detected without considering a priori data about SWE and LWP and quantifying their contributions. In particular, the following main questions are addressed: • To answer the above questions, we extract and isolate the contributions of snowfall scattering, snow cover scattering, and LWP emission in observed Tbs as a function of snowfall rate (sr), LWP, and SWE. To that end, we rely on multi-year coincident GMI [7] and DPR data [50], as well as ancillary information of SWE, LWP, and temperature. As previously noted, we mostly focus on the high-frequency channels at 89 and 166 GHz that are critical for snowfall retrieval.
Section 2 describes the products and data used for the analyses. Section 3 discusses the observed climatology of SWE, LWP, surface, and atmospheric temperature using their marginal and spatial distributions at different snowfall rates. By removing atmospheric effects, Section 4 quantifies the clear-sky microwave emissivity of the snow-covered surfaces as a function of SWE. In Section 5, we add the effects of LWP emissions, and finally, in Section 6, we complete the radiative budget by adding the snowfall component. Discussion and conclusions are presented in Sections 7 and 8.

Data and Products
The GPM core observatory, launched in 2014, carries the dual-frequency precipitation radar (DPR) and the GPM GMI, allowing for active and passive, colocated in space and time, observations [51,52]. Snowfall events and their radiometric signatures were extracted using the level-2 precipitation phase data from the GPM DPR product (2ADPR-V06) at the normal-sensitivity (NS) scan [53] and the calibrated GPM microwave imager (GMI) Tbs (1C-R GMI V05) [7]. Total precipitable water (TPW) is also from the GPM active microwave product (2ADPR-V06). The PMW measurements are from the GMI with 13 channels ranging from 10 to 183 GHz. In the DPR product, both Ka-and Ku-band (35.5 and 13.6 GHz, respectively) radar reflectivity values were used to estimate the precipitation rate to reduce the uncertainty of single band retrievals. The combination of active and passive observations provides a unique opportunity to understand unknown radiometric properties of solid precipitation [54] over snow-covered surfaces at frequencies above 100 GHz. The surface temperature (T skin ), 2-m temperature (T 2m ), the cloud LWP, vapor water paths (VWP), and ice water path (IWP) were all obtained from the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) [55]. The average air temperature, also derived from the DPR atmospheric state environmental ancillary information (2AD-PRENV), comes originally from the Japan Meteorological Agency (JMA) Global ANALysis (GANAL). This temperature is the average air temperature from 0 to 20 km at a 250-m resolution and will be called T air , hereafter in this paper. The extent of snow-covered surfaces was determined from the Interactive Multisensor Snow and Ice Mapping System (IMS) at 1 km spatial resolution [56], while SWE was obtained from the hourly MERRA-2 product [57] [M2T1NXLND] at a 0.625°× 0.5°spatial resolution. All data and variables used/calculated in this study, their units, and sources are listed in Table 1. In this paper, we only focus on events over dry snow cover, where both T skin and T 2m are below 0°C. Surface is considered as snow-covered when the IMS binary product indicates snow cover and the SWE from MERRA-2 reanalysis is greater than zero. We considered data from January 2015 to December 2020 over the Northern Hemisphere (NH) land surfaces. MERRA-2 data were linearly interpolated onto the scanning time of the GPM overpasses and re-gridded onto the DPR grid using the nearest neighbor interpolation. In the explained methodology, no information is lost at the expense of some redundancy in the coincident data. The binary snow-covered/no snow-covered IMS 1-km data was converted to 4-km binary snow-cover data using nearest neighbor (to find the IMS pixels in each 4 km by a 4-km 2ADPR grid) and the majority vote rule (>50%).

North Hemisphere Climatology of Snow Cover and Cloud Liquid Water
The radiometric signal of snowfall, snow cover, and cloud LWP are tightly intertwined at high microwave frequencies (≥89 GHz). The snowflakes and snow cover grains both scatter the upwelling surface emission (cooling effect), while the cloud LWP content adds to this emission (warming effect). An increase in SWE leads to an increased scattering volume of the snow cover and thus to a decrease of the surface emissivity and observed Tbs. An increased scattering of snowing clouds should further decrease the Tbs, especially for high snowfall rates. The snowfall scattering, however, becomes relatively less significant at larger SWE values because of a reduced background emission, which adds to the complexities of snowfall retrievals.

Marginal and Spatial Distribution of Snow Cover in Terms of Snow Water Equivalent
To unravel the effects of SWE and LWP on the quality of snowfall retrievals, we first need to investigate the frequency distribution of these variables. Figure 1a shows the positively skewed marginal probability distribution function (PDF) of SWE over the NH. The median and the mean are 41 and ∼74 kg m −2 , respectively. The histogram is color coded to show the mean T skin for each 25 kg m −2 SWE interval. SWE values between 30-300 kg m −2 show colder temperatures than the tail segments, with the lowest temperature occurring at about 75-100 kg m −2 . Figure 1a also shows that these SWE values with the lowest temperature on average occur at latitudes 54-57°N. It is calculated by averaging latitudes associated with the SWE data that fall into each of the SWE bins. To understand the reasons for the T skin behavior, we look into the maps of the SWE annual probability of occurrence and the corresponding mean T skin shown in Figure 1c-h. It is worth noting that here the mean refers to the entire NH in the histogram in each SWE bin, while it refers to the single grid box (1-deg resolution) in the maps.    (Figure 1c) are consistently displayed over land surfaces that are likely to receive snowfall (>30°N). The highest probability of occurrence is over the High Mountain Asia and southern Siberian Plateau, where T skin is visibly colder than the surrounding areas (<265 K, Figure 1f). However, the area with colder T skin (green area in Figure 1f) with a high occurrence frequency of 0-40 kg m −2 SWE values (red area in Figure 1c) is much smaller than the warmer T skin regions (orange and red areas in Figure 1f). Therefore, despite its higher probability of occurrence, the latter dominates the T skin mean values, and we see a relatively high T skin in the SWE PDF for low SWE values (0-40 kg m −2 ) in Figure 1a.
The  (Figure 1h). High values of SWE usually occur during late winter/early spring due to snow accumulation when the temperature begins to increase, justifying the slightly higher mean T skin compared to the 40-100 kg m −2 SWE interval. This higher mean temperature concurrent with larger SWE values (>100 kg m −2 ) over Pacific Coast Ranges and the Ural Mountains can weaken or mask the expected incremental increase in the scattering signal of snow cover in these regions.

Marginal and Spatial Distribution of Cloud Liquid Water
The parameterization of the cloud liquid water content over the snow-covered land surface requires an accurate characterization of the sub-grid distribution of thermodynamic variables such as cloud phase, cloud type, cloud vertical structure, precipitation occurrence, and geolocation. To this end, we investigate the marginal PDF and the spatial distribution of the multi-year LWP with respect to the snowfall occurrence as well as snowfall rate (sr) over dry snow-covered areas (defined as snow cover where both T skin and T 2m are below zero). The distribution of LWP as a function of T air and sr is presented in Figure 2. Similar to Figure 1, the mean refers to the entire NH in the histograms in each LWP bin, while it refers to the single grid box in the maps. Figure 2a-c shows the histogram of the LWP values corresponding to the orbit-level DPR observations using nearest-neighbor interpolation to match the 2ADPR-V06 resolution. The colors in the histograms shown in Figure 2a,b represent the multi-year average of T air for each LWP bin size of 14.2 g m −2 . For non-precipitating clouds (sr = 0, Figure 2a) and over dry snow cover-when the cold weather regime usually dominates-the PDF of LWP has a large positive skewness (γ = 1.91). This significantly skewed distribution indicates a high probability of occurrence of supercooled liquid clouds containing low LWP values. The mean and standard deviation are 36 g m −2 and 46 g m −2 , respectively. In the snowing atmosphere (sr > 0, Figure 2b), the PDF skewness decreases to γ = 0.08 and the mean increases to ∼140 g m −2 , indicating the existence of a larger amount of supercooled liquid droplets in clouds.
During snowfall events, the most frequent values of LWP are around 140-160 g m −2 , while the occurrence probability of extremely low LWP values (<20 g m −2 ) is below 0.0015 (Figure 2b). Larger values of LWP occur more frequently at warmer air temperatures as the moisture holding capacity of the atmosphere is higher and more ice water can turn into liquid [58]. However, the increasing rate of LWP with temperature is not constant. This is because different values of LWP occur at different vertical heights. The database developed by Kubota et al. [59] using the global cloud system resolving model for the GPM/DPR algorithm revealed that over land, the LWP increases with precipitation. Our analysis in Figure 2c confirms the evolution of the LWP probability distribution as a function of the snowfall rate. The PDFs become less skewed and wider as the mean moves from 80 g m −2 (sr = 0-0.5 mm h −1 , blue curve) to 200 g m −2 (sr = 4-8 mm h −1 , black curve). Furthermore, we observe that for sr > 2 mm h −1 , the LWP distribution does not change or shift noticeably. Figure 2. The effects of different overland snowfall rates (sr) on frequency and spatial distribution of cloud LWP over the NH: probability distribution functions of LWP for sr = 0 (a), sr > 0 (b) and 0 < sr < 8 mm h −1 (c); spatial distribution of the mean LWP (d-f), and the T air (g-i) for nonprecipitating (d,g), 0 < sr ≤ 0.5 mm h −1 (e,h) and sr ≥ 1 (f,i). The histograms are color coded to represent corresponding mean T air temperatures. Note that the color represents the multi-year average of T air for each LWP bin (in the top row), the multi-year gridded average of LWP (second row), and the multi-year gridded average of T air (bottom row).
To understand the spatial variations of LWP from non-snowing to snowing atmospheres, we stratified our dataset based on different snowfall intensities. Figure 2 shows the spatial pattern of mean LWP (d-f) and T air (g-i) for three different scenarios: sr = 0 (no snowfall), 0 < sr ≤ 0.5 mm h −1 (about 40 percentile), and sr ≥ 1 mm h −1 (95 percentile)over dry snow cover. The average LWP and its corresponding T air increase with the increase of sr. This is because the predominant snowfall-temperature relationship is positive at mid-to high-latitudinal regions during the cold-weather regime in winter [60]. The snowfall warming feedback significantly increases the T air over North America, Siberia, and East Asia where colder temperatures were observed in the absence of precipitation, while temperatures at lower latitudes remain almost unchanged (Figure 2g-i). Figure 2 also shows that the rate of increase in both LWP and T air with snowfall rate is not uniform over the entire NH. The increase of T air specifically is observed over the Siberian plateau and northern Canada. The increase of LWP over these regions concurrently occurs with an increase of atmospheric temperature, both of which increase the emissivity and could significantly mask the Tb response to the increase in snowfall scattering.

Radiometric Effects of Snow Cover and Cloud Liquid Water on GPM Brightness Temperature
To understand the effects of the observed complex climatology of snow cover, temperature, LWP, and morphology of snow cover grain on GPM brightness temperatures, the averages of high-frequency Tbs are calculated as a function of SWE at different LWP values. Figure 3 shows the average NH Tbs as a function of SWE and cloud LWP in the absence of precipitation. For brevity, we focus on the vertical polarization (V-pol) channels as the response pattern is similar in the horizontally polarized (H-pol) ones (not shown here). We would expect a monotone decrease in the snow cover emissivity pattern as the SWE increases for frequencies below 100 GHz [23,24,48]. However, results show that Tb values do not monotonically decrease and there is a clear inversion in the Tb spectra. As shown in Figure 3, Tb values approach a minimum at around 70-100 kg m −2 and begin to increase for SWE > 120 kg m −2 . This anomalous spectra was previously observed in SSM/I frequency channels at 19, 37, and 85 GHz [30] and explained as a radiometric response of the snowpack to its crystalline structure changes when the mean grain size increases [27]. Among the GMI channels, both the 89 and 166 GHz ones are almost equally sensitive to the increase of SWE, while the 183 GHz is the least sensitive, because of the atmospheric water vapor emission. On average, for a 1 kg m −2 increase in SWE, the Tb decreases about 0.5 K at 89 and 166 GHz, 0.2 K at 183 ± 3 GHz, and 0.4 K at 183 ± 7 GHz. It is worth noting that the observed anomaly follows the climatology of T skin for different SWE values (Figure 3f). The snow cover metamorphic changes on grain size and the snow cover temperature climatology are not independent; therefore, the observed scattering signal is a response to both seasonal and spatial variations of snow cover emissivity and T skin , as shown in Figure 1.
Since the crystalline structure of the snow pack is a function of its metamorphic changes, the observed anomaly shows a seasonal dependence. The Tb seasonal dependence here is calculated by averaging the time stamp of Tb values at each SWE bin. Figure 3 shows that the maximum snow cover scattering occurs during the early winter. As the T skin and the snow mean grain size increase toward the late spring, the scattering signal begins to decay even though the SWE continues to increase. It is important to note that this seasonal pattern impacts the quality of snowfall microwave retrievals. In particular, since strong snow cover scattering can weaken the snowfall signal, we expect larger uncertainties on PMW snowfall retrievals on early winter when the snow cover is fresh and the SWE is less than 100 kg m −2 .
The analysis also shows warmer Tbs for increasing frequencies which is in contradiction with the known surface emissivity spectra of snow cover [23]. Within the analyzed range of SWE (0-400 kg m −2 ), the 166 GHz channel is more than 5 K warmer than the 89 GHz one, while the overall expectation is that the high-frequency channels must be colder due to stronger snow cover scattering [23,24,48]. This inverted spectrum was observed in lower frequency channels of the SSM/I sensor and is known to be due to the formation of a dense layer of snow crust [30]. Hewison and English [61] also argue that this phenomenon could be due to the mean behavior of the snow cover temperature profile, which is often colder at the bottom layers. Finally, an increase of cloud LWP can completely mask the scattering effects of snow cover. This masking effect is shown in Figure 3 as we see the curves are flatter for higher LWP (moving from blue to red lines). This is expected because of the climatology of LWP and T air observed in Figure 2, e.g., larger values of LWP occur at warmer T air , on average. Figure 3a-d shows that for LWP > 150 g m −2 , there is almost no response to changes on SWE. Moreover, the SWE value associated with the maximum scattering (Tb minimum) increases as the cloud LWP increases. The reason could be while the maximum snow cover scattering occurs in early winter (Oct-Dec) over dry snow, the emission of cloud liquid water is stronger from mid to late winter (Jan-Mar) [62]. This is attributed to the fact that the temperature begins to gradually increase from mid to late winter. This masking effect is also frequency-dependent. In particular, when the scattering signal of snow cover is maximum, the Tb 89 GHz increases about 5 K more than the Tb 166 GHz in response to an increase of 200 g m −2 in cloud LWP. We also observe that the difference between the channel spectra shrinks for high (200 kg m −2 ) and low (2.5 kg m −2 ) values of LWP (Figure 3e)-revealing the nonlinear effects of LWP on snow cover high-frequency signatures. Now that we have discussed and showed the complex radiometric effects of SWE and LWP climatology and snow-cover grain morphology on the GPM brightness temperatures, the question is how we can isolate these correlated effects on the Tb signal of snowfall in order to improve snowfall passive retrieval from GPM frequency channels. We approach this by calculating emissivity in the following sections.

The Snow Cover Emissivity under Clear Sky
The microwave signal reaching the top of the atmosphere is the combination of the surface signal and the atmospheric contributions. If the atmospheric contribution becomes small enough to be ignored, we can easily calculate the surface emissivity (ε s ) as follows, as also used by [63,64] and evaluated by [65]: where Tb s is the brightness temperature observed at the top of the atmosphere in clear-sky conditions, and thus, it accounts only for snow-cover land. To minimize the atmospheric contribution, we only consider clear-sky observations with zero LWP, zero IWP, and zero precipitation. The atmospheric WVP is still present and thus might affect ε s in the 183.3 ± 3 and 183.3 ± 7 GHz channels [41] but can be ignored for other frequencies. Figure 4 shows the emissivity of snow-covered surfaces calculated for clear-sky as a function (Equation (1) Figure A1 (Appendix A), are very similar to what previous emissivity studies found [23,66]. Figure 4a,b shows that for SWE values smaller than ∼100 kg m −2 , the sensitivity of the surface emissivity to the SWE increase is larger at 89 GHz compared to 166 GHz. This trend is reversed for SWE ≥ 100 kg m −2 , with the 89 GHz emissivity reaching a plateau at around 0.83 for H-Pol and 0.86 for V-Pol. This indicates that the emissivity at channel 166 GHz provides a better response to a full dynamic range of SWE with a consistent ε s decrease, while at 89 GHz, it saturates at SWE ∼100 kg m −2 thus almost becoming blind to any further increase in SWE. We expect that this blindness of channel 89 GHz to SWE mostly occurs over the Pacific Coast Ranges of northern British Columbia, North East Canada, Ural Mountains, Kjolen Mountains of Norway, and Verkhoyansk range in Russian Far East (Figure 1e) during Jan-May ( Figure 3).
The variation of ε s at water vapor channels (183.3 ± 3 and 183.3 ± 7 GHz) is very small (∼0.015) and strongly anti-correlated with the variations of the atmospheric WVP (Figure 4f). This is obviously due to the effects of the WVP in the atmosphere masking the surface contribution. Therefore, hereafter in this paper, we only focus on the high-frequency window channels (89 V, H and 166 V, H GHz).
Our key observations on emissivity change rates in response to the increase in snow depth are summarized over the dry snow-covered surfaces after removing atmospheric contributions: (1) for SWE less than 10 kg m −2 , the decreasing rate of emissivity at 89 GHz exceeds that at 166 GHz and on the contrary, (2) for SWE greater than 100 kg m −2 , the snow cover emissivity decays at a much higher rate at 166 GHz compared to 89 GHz. Although the 89 GHz emissivity of the surface covered with dry snowpack reaches a plateau for SWE ≥ 100 kg m −2 , it remains significantly lower than the non-plateaued emissivity at 166 GHz. This low emissivity can better capture the emission signal of snowfall, which we analyze in the next two sections. , and 183 ± 7 (e)), the emissivity differences of V-pol and H-pol at 89 and 166 GHz frequencies (c), and the corresponding WVP (f). Note that the total cloud liquid water path, integrated total precipitable water, and ice water path in the column of atmosphere and the near-surface precipitation are all zeros based on the data from both MERRA-2 model simulations and the DPR ancillary database (2ADPR-ENV).

Effects of Cloud Liquid Water Emission on Snow Cover Emissivity
Tb values at high-frequency channels may increase during snowfall events if clouds contain some liquid or supercooled water, which could completely mute the scattering contribution of snowfall [6,22,38,67,68]. Therefore, it is crucial to quantify the mixing effects of SWE scattering and thermal emission of LWP on the snowfall signal. We calculate the emissivity for different combinations of SWE and LWP in the absence of precipitation. We compute the emissivity due to the LWP emission, assuming that the observed Tb (Tb obs ) is the sum of the clear-sky Tb (Tb s = ε s ·T skin from Equation (1)) and the Tb due to the atmospheric LWP (Tb lwp ) as: Tb obs = Tb s + Tb lwp Tb obs = ε s T skin + ε lwp T air (2) where ε lwp is the emissivity due to the cloud liquid water, ε s is the snow-cover emissivity under clear sky, and T air is the air temperature. Therefore, the emissivity due to the LWP emission can be calculated as follows, according to Equation (2): The calculated ε lwp by Equation (3) are presented in Figure 5. Here, ε s is calculated for each SWE value using the regressed relationships shown in Figure 4  The increase of LWP (from bottom to top curves) increases the atmospheric emissivity, and this variation can be significantly better captured at 89 than 166 GHz ( Figure 5). This is because the 89 GHz clear sky emissivity (Figure 4a,b) is lower than the 166 GHz one, which allows the 89 GHz channel to have a more pronounced response when adding the emission contribution of LWP over SWE compared to channel 166 GHz. However, Figure 5a shows that this contribution progressively decreases and becomes very small (∼0.005) for LWP ≥ 200 g m −2 . These large values of LWP as observed in Figure 2d occur mostly over the Appalachian Mountains, the west coast of British Columbia, Alaska, and Northern Europe. On the contrary, at channel 166 GHz, the LWP emissivity contribution shows a noticeable increase (∼0.03) only for LWP values up to about 100 g m −2 , while for higher values, LWP needs to increase significantly to show a similar emissivity response (Figure 5b). This observation is very critical for snowfall detection since the LWP, as observed in Figure 2c, largely increases during snowfall events. Moreover, larger emission contributions from clouds could mask the variations of surface emissivity caused by different snow cover scattering, particularly at the 166 GHz channel. We see both channels respond better to the scattering of snow cover for SWE values less than 100 kg m −2 . This is obvious in Figure 5 as lines are steeper for LWP when SWE is less than ∼100 kg m −2 . Similar findings for WVP were also revealed by [18] about the water vapor increase in the column of the atmosphere that could mask the variations of radiometric response to surface temperature at frequencies ≥ 89 GHz.

The Interactions of Snowfall Scattering with Snow Cover and Cloud Liquid Water
Now that the SWE and LWP emissivities are quantified, the snowfall scattering signal can be isolated as follows: Tb obs = Tb s + Tb a Tb obs = ε s T skin + ε a T air ε a = ε sr + ε lwp (4) where Tb a and ε a are the atmospheric brightness temperature and atmospheric emissivity, respectively. Furthermore, recall that Tb s from Equation (1) is the Tb at the top of the atmosphere at clear-sky, which only accounts for snow-cover emissivity.
The calculated ε a from Equation (4) is shown in Figure 6. It is worth noting that when sr = 0, ε a and ε lwp are equal. The difference between emissivities in Figure 6 for sr = 0 and sr > 0 shows the snowfall scattering contribution reducing the ε a as all of the curves have almost equal amount of liquid water content (150-200 g m −2 ). This LWP interval is selected based on its high probability of occurrence shown in Figure 2 to ensure an adequate sample size in each SWE interval (20 kg m −2 ) and snowfall rate interval of 0 < sr ≤ 0.5, 0.5 < sr ≤ 1, and 1 < sr ≤ 4 mm h −1 . The reason we did not go beyond sr = 4 mm h −1 is that the sample size significantly reduced due to the lack of enough heavy snowfall events in each LWP and SWE interval. . Note that the dark gray curve is associated with sr = 0 mm h −1 , which is the exact same curve (solid dark gray with yellow markers) in Figure 5. It is shown here again for making easier comparisons of ε a when sr > 0 and ε lwp .
The results in Figure 6 show that for SWE values close to zero and the same amount of LWP, the snowfall signal with intensity 0-0.5 mm h −1 reduces the ε a at 89 GHz by about 27% (from 0.11, grey line to 0.08, red line- Figure 6a) and at 166 GHz by ∼38% (from 0.09, grey curve to 0.055, red curve- Figure 6b). At 89 GHz, this decrease in emissivity due to snowfall scattering decreases with SWE, with a gap of about 0.015 (∼18%) for SWE ∼200 kg m −2 and less than 0.005 (∼3%) for SWE greater than 400 kg m −2 ( Figure 6a). As illustrated in Figure 4, for the clear sky emissivity at the 89 GHz channel when the snow depth increases on the ground, the land surface emissivity becomes relatively small (0.83) and reaches a plateau. Therefore, the small snowfall scattering cannot reduce the emissivity any further. However, the snowfall with sr 1 < sr ≤ 4 mm h −1 can further decrease the emissivity by about 0.018 from 0.078 to 0.06 for SWE values of <40 kg m −2 and from 0.03 to about 0.01 at SWE values larger than 400 kg m −2 (the difference between red and blue curves, Figure 6a).
For channel 166 GHz, Figure 6b shows that the snowfall scattering contribution (refers to the emissivity decreases due to snowfall scattering) is almost the same as that of 89 GHz when the intensity is small (27%, sr ≤ 0.5 mm h −1 ). However, channel 166 GHz responds more strongly to the increase in snowfall rate showing clearly separated curves for different sr intervals. The snowfall scattering contribution when 1 < sr ≤ 4 mm h −1 can further decrease the emissivity by about 0.06 (difference between the red and blue curves ∼96%) at channel 166 GHz compared with only 0.018 at channel 89 GHz (∼23%).

Discussion
By integrating the findings from Figures 5 and 6, we conclude that for small snowfall intensities (≤0.5 mm h −1 ) over shallow snow cover (SWE < 100 kg m −2 ), LWP values even as small as 0-50 g m −2 could increase the emissivity by 0.02-0.05 and thus, could completely mask the scattering contribution of snowfall (∼0.05 as a difference between ε lwp and ε a ) at channel 89 GHz. Larger snowfall rates might only be captured at small LWP values. The 166 GHz channel could still capture some snowfall scattering at SWE values less than 100 kg m −2 , defeating the emissivity increase of small liquid water content. However, the snowfall scattering contribution is also masked at this channel when the LWP becomes larger than 100-150 g m −2 , which increases the emissivity by almost more than about 0.07.
Over deeper snow cover with SWE larger than 100 kg m −2 , the snow cover scattering contribution becomes very significant and thus, alleviates the contribution from the small intensities of snowfall scattering, making them no longer distinguishable from the background emissivity at both 89 and 166 GHz channels. Larger snowfall intensities (1-4 mm h −1 ) still decrease the background surface emissivity even over deep snow cover with SWE ≥ 100 kg m −2 at the 89 GHz channel for LWP values up to 50-100 g m −2 but are masked for larger LWP values. This is because the decrease in emissivity due to the increased snowfall rate (Figure 6a) is lower than the increase in emissivity because of the increased emission of LWP (Figure 5a). At 166 GHz, both SWE scattering and LWP emission over deep snow-covered surfaces are relatively smaller than those at 89 GHz. Therefore, the 166 GHz channel is more capable of capturing the snowfall scattering signal with large intensities (sr > 1 mm h −1 ) for LWP values up to 100-150 g m −2 .
For larger LWP values (∼>150 g m −2 ) over regions with SWE ≥ 100 kg m −2 , it is revealed that the snowfall signal could be captured by its emission signature-LWP increase during snowfall events-instead of its scattering signature at channel 89 GHz. This is because the surface emissivity at channel 89 GHz reaches a plateau value and remains unchanged to further increases in SWE. Therefore, over regions with deep snow cover, the 89 GHz channel could capture this emission signal of the snowy clouds for snowfall detection. The challenging land-atmospheric situations for snowfall retrieval particularly occur if there is not any scattering nor emission signal, which is when LWP > 100-150 g m −2 and SWE < 100 kg m −2 .
In the present article, the relations between Tbs, snowfall, cloud liquid water with emphasis on snow-covered regions were established as multi-year averages from four years of colocated GPM and MERRA-2 reanalysis data. The immediate goal here was to investigate the challenging zones of snowfall retrieval over the snow cover at high-frequency GPM channels regarding the confounding effects of these atmospheric constituents over the snow-covered regions. However, at the daily or sub-daily time scales, large variability around these multi-year averages is expected. When seeking to establish relationships at these finer time scales, the temporal variability needs to be either handled as a stochastic process or analyzed in relation to additional physical parameters, which have not been taken into account in the present study (e.g., snow density and wetness, particle shape and size distribution). Additionally, at fine time scales, significant variability is expected to be related to errors and inaccuracies in the measured/estimated variables in both satellite observations and reanalysis products.
For instantaneous precipitation retrievals, one should consider the use of the dynamic surface emissivity database developed by [69] and evaluated in GPM retrieval by [70], which implements the optimal estimation method with a forward model error covariance matrix. Future research would benefit from using this dynamic land surface emissivity database for further investigation of the findings of the present paper and for establishing relationships between LWP, SWE, and snowfall at finer temporal scales.
The fact that GMI and the DPR are flying onboard the same platform allows for a high number of colocated observations at all latitudes; however, we acknowledge potential errors and inaccuracies in DPR measurements regarding light precipitation intensities. An additional investigation, if it does not require high space-time coverage, should consider measurements from CloudSat Cloud Profiling Radar (CPR), which is known to more accurately capture low-intensity snowfall [71,72].
In this study, we only used total precipitable water to screen the clear-sky Tbs. This can add some uncertainties regarding the calculated emissivities of LWP and snowfall. Future research needs to investigate the effects of total precipitable water on the radiometric signal during the snowfall events [73].

Conclusions
We analyzed snow cover effects on the snowfall radiometric signal under different snow water equivalent and liquid water content of clouds, by considering the effect of surface temperature variations, as well. To isolate the surface temperature thermal emissivity, we first calculated the clear-sky emissivity. Then, we computed the contribution of scattering and emission from the other mentioned land-atmospheric constituents directly by comparing their emissivity with clear-sky snow cover emissivity. We found that the channel 166 GHz could better capture the scattering signature of light snowfall events because it responds less strongly to the increase of the cloud liquid water than the 89 GHz channel. Larger snowfall events could be captured better at 89 GHz when both LWP and SWE are small, while 166 GHz becomes more advantageous at capturing this scattering when LWP increases up to about 100-150 g m −2 . Over deeper snow-cover regions, and particularly larger LWP values (≥100-150 g m −2 ), the scattering of snowfall, even with large intensity, is masked by the comparable scattering contribution from the large accumulation of snow cover and the emission from liquid water at both 89 and 166 GHz channels. At this land-atmospheric condition, the snowfall dominant signature becomes its emission that can be distinguished from the very low plateaued emissivity of the surface at channel 89 GHz. We believe that a quantitative climatological assessment such as this presented herein can provide useful information for improving passive microwave retrieval of snowfall and also serve as diagnostics for interpreting the bias and uncertainty of current products.  Data Availability Statement: The GPM data are provided courtesy of the NASA Precipitation Processing System at the Goddard Space Flight Center (https://pmm.nasa.gov/data-access/, access date: 24 June 2021). The MERRA-2 data are from the Goddard Earth Sciences and Information Service Center (https://disc.sci.gsfc.nasa.gov/mdisc/, access date: 24 June 2021). The IMS snow cover product is available at United States National Ice Center (ftp://sidads.colorado.edu/pub/ DATASETS/NOAA/G02156/, access date: 24 June 2021). All colocated GMI, SWE, LWP, WVP, TPW, sr from DPR, T skin , T air , and T 2m at DPR resolution and calculated emissivity data ε s , ε LWT , and ε a are available upon request.

Acknowledgments:
The authors would like to thank F. Joseph Turk from the Jet Propulsion Laboratory, Gail Skofronick-Jackson from NASA HQ and Ardeshir Ebtehaj from St. Anthony Falls laboratory at the University of Minnesota for their valuable feedback at the early stage of this research.  Figure A1. Clear-sky emissivity with snow water equivalent at GPM low-frequency channels (10.6, 18.7, and 36.5 GHz V, H) and the emissivity differences of V-pol and H-pol. Note that the total cloud liquid water path, integrated total precipitable water, and ice water path in the column of the atmosphere, and the near-surface precipitation are all zeros based on the data from both MERRA-2 model simulations and the DPR ancillary database (2ADPR-ENV).