Abstract

A physically based method is developed to estimate the microphysical structure of the melting layer in stratiform rain using airborne observations by a dual-frequency radar and a 10.7-GHz radiometer. The method employs a nonlinear optimal estimation approach to find two parameters of the gamma drop size distribution (DSD) at each radar range gate from the Ku/Ka-band reflectivities. The DSD profile is used to determine the atmospheric absorption/extinction profile, which enables the surface contribution to the measured brightness temperature to be estimated. The surface wind speed is estimated from the surface emissivity by inverting the forward model, which relates the two. Retrievals in stratiform precipitation require a model to describe the thermodynamic and electromagnetic properties of melting hydrometeors. The melting layer can contribute a majority of the total atmospheric absorption, making it a key component for accurate retrievals in stratiform rain. Several melting layer models were evaluated based on their fit to the dual-frequency reflectivity measurements in the melting layer. A candidate model is selected and tuned to match the radar measurements. The melting layer model is then incorporated into the full forward model for the brightness temperature observed by the radiometer. The surface wind speed assumed in the forward model is forced by the radiometer observations. If the actual surface wind speed is known, this approach provides a powerful constraint on the possible melting layer model. A case study is presented from an airborne campaign over areas of precipitation off the coast of Vancouver Island, British Columbia, Canada. The estimated wind speeds are found to be uncorrelated with the reflectivity and their average value is within 1 m s−1 of that retrieved in a clear area adjacent to the rain.

1. Introduction

Knowing the surface wind field under areas of precipitation could lead to better prediction of severe weather phenomena and how they will affect the human population. For example, the accuracy in the prediction of hurricane storm track and intensity is extremely dependent on the accuracy of the wind measurements within the hurricane. In recent years, hurricane prediction models have vastly improved in their ability to predict storm track, but predictions of storm intensity have proved to be much less reliable and are an area of ongoing research. Combining accurate estimates of the surface wind field with vertical profiles of the precipitation water content, and associated latent heat release, will be invaluable for use in hurricane forecast models. Remote sensing instruments are seen as a valuable tool for making these measurements because they provide extensive spatial coverage and make measurements where it is impractical or physically impossible to have an in situ instrument.

One impediment to accurate remote sensing retrievals of surface wind and precipitation is the inadequacy of current models for the melting layer present in stratiform rain. The objective of this study is to develop a physically based method for developing such a melting layer model, and for validating it that is constrained by microwave radiometer and radar observations in areas of convective and stratiform precipitation over the ocean. Investigations of the melting layer have increased over the past decade due to the emergence of physically based precipitation retrieval algorithms. It was shown that neglect of the excess attenuation in the melting layer could cause errors as high as 100% in passive microwave retrievals of precipitation (Bauer et al. 1999). It is shown in this paper that neglect of the melting layer would lead to errors much greater than this for a retrieval algorithm that is attempting to retrieve surface parameters in rain, such as wind speed.

Our development of a melting layer model is imbedded in a physically based retrieval algorithm, in which parameterized characteristics of the melting layer are adjusted along with rainfall and surface roughness parameters in order to best fit the microwave observations. Previous combined radar–radiometer algorithms have only attempted to retrieve rainfall parameters and assumed a fixed surface wind speed. Current methods to measure the surface wind speed under rain only work at wind speeds greater than 10 m s−1 and are not able to retrieve the vertical profile of precipitation, only providing a path-averaged estimate (Uhlhorn and Black 2003). A new method is developed in this study that retrieves both the near surface wind speed and the vertical precipitation content from an airborne 13.4- and 35.6-GHz radar and a 10.7-GHz horizontally polarized radiometer.

The algorithm developed in this study solves for two parameters of the gamma drop size distribution at each radar range gate by employing an iterative inversion technique using Mie scattering theory. Retrieving the drop size distribution directly will eliminate most errors associated with Z–R or k–R relations in which a fixed size distribution (DSD) is required, since the DSD is known to change from storm to storm and over the course of a storm’s evolution. The parameterized DSDs are then used to determine the atmospheric extinction and absorption present in the measured brightness temperatures in order to isolate the surface emissivity, which has a monotonic dependence on the near surface wind speed.

In light stratiform rain, where the freezing level (FL) is close to the surface, the most important component of the retrieval algorithm is the model that describes the microphysical properties and dielectric constant of hydrometeors in the melting layer, and thus predicts their absorption and scattering. Previous models for the melting layer exist, but they tend to predict vastly different amounts of scattering and absorption. Currently, there is no general agreement on which model to use. Models that over or underestimate the absorption in the melting layer will, respectively, under or overestimate the retrieved wind speed in stratiform rain. Different melting layer models are evaluated in this study based on how well they are able to predict the reflectivity profile measured at two radar frequencies. A candidate model is selected and tuned to fit the radar observations. This updated model is then used in the retrieval algorithm and is shown to estimate physically reasonable wind speeds in stratiform rain that contains varying intensities of peak reflectivity in the melting layer. The high sensitivity of the retrieved wind speed to the melting layer model presents a novel technique to constrain the melting layer. An analysis is performed to determine with what fidelity one could separate various melting layer models if the surface wind field were known to a given accuracy. This provides a third sensitive constraint on the melting layer, in addition the dual-frequency reflectivity profiles.

A case study is presented that demonstrates the retrievals in stratiform and convective rain. The data was acquired from a field campaign in which Jet Propulsion Laboratory’s (JPL) Second Generation Precipitation Radar (PR-2) Ku–Ka-band radar and the University of Michigan/Goddard Space Flight Center Lightweight Rainfall Radiometer X-band radiometer (LRR-X) were flown over areas of precipitation off the coast of Vancouver Island, British Columbia. A thorough sensitivity analysis is developed to estimate the uncertainty in the retrieved products and to suggest future instrument configurations to improve the retrievals.

2. Instrumentation

a. Lightweight Rainfall Radiometer

The experimental LRR-X was designed and built at the University of Michigan and the National Aeronautics and Space Administration (NASA) Goddard Space Flight Center (Ruf et al. 2002; Ruf and Principe 2003). LRR operates at 10.7 GHz with horizontal polarization. LRR is a synthetic thinned aperture radiometer that has 14 slotted waveguide antennas on a 1-m2 aperture, optimally placed to measure the complex cross correlation between pairs at 68 λ/2 spacings. By sampling the observed interference pattern, the brightness temperature image is able to be reconstructed at very high resolution cross track. The antennas have a 2.1° angular resolution along track. At the DC-8 altitude (∼10.5 km), LRR has a pixel resolution of 381 × 466 m at nadir, increasing to 1079 × 629 m at 45° cross track. The NEΔT of the system is 0.3 K.

The calibration of LRR consists of two parts. First, the raw voltage outputs from the correlators are calibrated using a reference load and noise diode source. This produces 68 visibility samples of the scene. The visibility samples are converted to the brightness temperature image using a discrete inverse Fourier transform. A final end-to-end TB calibration of each pixel in the image is necessary to account for antenna sidelobes and hardware loss between the antennas and the receiver. The LRR brightness temperatures are calibrated to modeled TBs over clear-sky areas encountered immediately before and after the precipitation overflights. Linear interpolation of the clear-sky calibration parameters are used to calibrate the data in between the two calibration points to account for system drifts. The ocean brightness model used for calibration is the same as that used in the retrievals. A coincident radiosonde sounding was used to determine the atmospheric absorption profile. The sea surface temperature and wind speed in the clear-sky areas were determined from a nearby buoy.

b. Second Generation Precipitation Radar

The PR-2 was built and is operated by the Jet Propulsion Laboratory in Pasadena, California (Im et al. 2002). The PR-2 measures the reflectivity at 13.4 and 35.6 GHz and mechanically scans to ±25° cross track. It also measures the Doppler radial velocity at 13.4 GHz and the linear depolarization ratio at 13.4 GHz. The vertical resolution is 37 m, and the horizontal resolution at the surface is 800 m from the DC-8 altitude. It has a 0.4-m offset reflector antenna, which is underilluminated at 35.6 GHz, producing matched beams at both frequencies. The minimum detectable reflectivity at both frequencies is approximately −15 dBZ at a range of 2 km from the radar and approximately 0 dBZ at 10 km. Contamination is observed near the surface due to the pulse compression sidelobes. The surface contamination begins approximately 600 m above the surface and must be accounted for, especially in light rain. When the reflectivity of the rain near the surface is higher than about 30 dBZ, the contamination is only evident beginning around 300 m above the surface.

c. Pacific field campaign

LRR and PR-2 were integrated on a NASA DC-8 aircraft and several science flights were made in June 2003. A midlatitude cyclone that formed in the Gulf of Alaska provided an opportunity to conduct precipitation overflights on 13 June 2003. Regions of stratiform and convective rain were observed. The freezing level was at approximately 1500 m and a strong bright band was observed in the stratiform regions. This mission was a proof of concept for the LRR technology, so no dedicated ground truth was acquired directly below the underflight regions. Although, surface measurements of the continuous and gust wind speed were taken by two buoys in the area of the overflights. The buoys are both maintained by Environment Canada. The first buoy is the Dellwood station at 50.86°N, 129.91°W and the second buoy is the La Perouse Bank station at 48.83°N, 126.00°W. This places the Dellwood station about 150 km north of the center of the overflight area and the La Perouse Bank station about 75 km south and 220 km east of the center of the observation area. The Dellwood station observed continuous winds near 11 m s−1 with gusts of 14 m s−1 during the observation time period. The La Perouse Bank station observed winds near 6 m s−1 with gusts of 7 m s−1. The sea surface temperature measured by the buoys at the time of the overflights was 284.1 K at the Dellwood station and 285.7 K at the La Perouse Bank station. A constant SST of 284.1 K was assumed for all analysis.

3. Algorithm to estimate DSD

The equivalent radar reflectivity, Ze(f) (mm6 m−3), of a volume of scatterers is related to the drop size distribution, N(D) (m−3), by

 
formula

where τf(0, R) (Np km−1) is the one-way path integrated extinction between the radar and R due to precipitation, cloud water/ice and atmospheric gases, K is a dielectric factor, N(D) is the number of drops of diameter D, and ξb is the Mie backscattering efficiency, which is a function of the index of refraction, n, of the particle, the particle diameter and the frequency of the incident wave.

This study assumes a normalized gamma DSD, which has been used by previous investigators (Willis 1984; Viltard et al. 2000). The normalized gamma DSD is similar to the traditional gamma DSD (Ulbrich 1983) with a normalization factor introduced

 
formula

where N*0 is the normalized offset term (m−4) and

 
formula

where Γ is the gamma function. The distribution of rain drops with size in a certain volume using the gamma DSD can be described by three parameters, N*0, D0, and μ. The rain rate and precipitation water content are proportional to the third moment of the DSD and can be computed given estimates of the DSD parameters (Rogers and Yau 1996 208–209). The expression from Atlas et al. (1973) is used for the terminal velocity of a given diameter rain drop when calculating the rain rate from the DSD.

The dual-frequency backscattering measurements can be used to estimate two parameters of the normalized gamma DSD at each radar range gate. The inversion is accomplished by an iterative nonlinear optimal estimation method (Rodgers 2000, 81–90). Values of N*0 and D0 are determined which minimize the RMS difference between the measured and modeled reflectivity at each range gate Rj using the estimator

 
formula

where x = [N*o, D0]T, y = [Z( f1), Z( f2)]T, the forward model, F(x) is given by Eq. (1) and the Jacobian, 𝗝, is calculated numerically. In this study, μ is set to a value of 0.99 that is constant with height and location. To better condition 𝗝, the value of N*0 is scaled by 10−6 for the inversion to have a similar magnitude as D0. Convergence is reached when (𝘆 − F(𝘅))T(𝘆 − F(𝘅)) < E, where E is an acceptable residual error. The measurement error covariance matrix, 𝗦E, is most simply represented as a diagonal matrix, with the variance of each measurement on the diagonal. The measurement error is assumed to be 0.5 dB of random additive Gaussian noise. Here 𝗦a is the a priori covariance matrix and represents the natural variability of the a priori state vector.

Although the inversion attempts to estimate two unknowns given two measurements, the backscatter measurements are not completely independent. A priori information about the DSD parameters is used to increase the stability of the inversion procedure and to prevent the retrieved parameters from diverging to nonphysical values. This is especially important when the rain rates become large and the attenuation is high. The a priori information used is xa = {N*0a, D0a}, which is the average a priori value of the two DSD parameters, and 𝗦a, which is a 2 × 2 covariance matrix of the a priori distribution of the two parameters. These values are set using the mean and variance of a uniform distribution of each parameter between physically allowable bounds. The maximum allowable bounds for D0 are 0.05 to 0.22 cm for stratiform rain and 0.05 to 0.26 cm for convective rain. The bounds for N*0 are also set based on whether the profile is convective or stratiform. The a priori bounds for stratiform rain are 1 × 106 to 8 × 106 m−4 and the bounds for convective rain are 13 × 106 to 25 × 106 m−4. For a given liquid water content, the stratiform DSD will have a higher concentration of large drops compared to the convective DSD and likewise the convective DSD will have a larger concentration of small drops compared to the stratiform DSD. These bounds bracket the offset parameters used by many investigators in radar precipitation retrieval algorithms. Viltard et al. (2000) sets N*0 to 17 × 106 m−4 for convective rain and 4.8 × 106 m−4 for stratiform rain. Dou et al. (1999) derived values from airborne microphysical probe measurements and found N*0 to be 19.2 × 106 m−4 for convective rain and 3.9 × 106 m−4 for stratiform rain. Grecu and Anagnostou (2002) use N*0 equal to 20 × 106 m−4 for convective rain and 2.2 × 106 m−4 for stratiform rain in their rain profiling algorithm. The bounds on N*0 in the snow layer are the same as those for stratiform rain, which was observed to be valid for snow particle distributions by Sekhon and Srivastava (1970).

The two-way attenuation at each range gate is estimated using the Hitschfeld–Bordan (1954) approach. In this approach, the attenuation from precipitation is assumed to be zero at the closest range gate to the radar with a measurable returned power, defined as range gate R1. The solution of the drop size distribution at R1 is used to determine the extinction coefficient due to precipitation, which is summed with the extinction due to gases and cloud liquid to determine the integrated attenuation at the next range gate, R2. This process is repeated from the first to the final range gate.

The Mie backscattering efficiency is computed based on whether the backscattering volume is in the rain, snow or melting layer. The formulation used in this study to calculate the Mie efficiencies is outlined in Ulaby et al. (1981). In the rain layer, spherical drops are assumed and the dielectric constant of freshwater is used for the rain drop (Ulaby et al. 1986), which is a function of the temperature of the drop and the frequency of the incident wave. For computational efficiency, the extinction, absorption and backscattering efficiencies are calculated at drop radii of 0.01 to 10 mm in 0.01-mm steps and at temperatures of 273 to 300 K in 1-K steps. These values are placed in lookup table for the computation of the volume extinction, absorption, and backscattering coefficient for a given number distribution of drop sizes.

The Mie efficiencies in the melting layer are tabulated as a function of height below the freezing level in 10-m increments based on the melted mass fraction of each particle predicted by a thermodynamic/microphysical model and the dielectric of each melting particle predicted by a dielectric mixing model. Details of these models are given in section 5. The Mie efficiencies for snow particles are taken as the special case of a melting particle with a melted mass fraction of zero. This means that the same dielectric mixing formula and snow density model will be used to determine the Mie efficiencies in the snow layer and in the melting layer.

A brightband detection algorithm is used to determine if a melting layer exists for a given reflectivity profile and, if so, to determine the height of the top and bottom of the melting layer. The algorithm works by exploiting the characteristic shape of a reflectivity profile containing a bright band. The algorithm is adapted from one described in Klaassen (1988). If no melting layer is found to exist, the profile is labeled as convective. For a stratiform profile, each range bin is designated as belonging to either the snow layer, the melting layer, or the rain layer. For a convective profile, each range bin belongs to either the snow layer or the rain layer, and the transition occurs at the freezing level determined from the algorithm.

4. Retrieval of wind speed

The brightness temperature of a scene observed by a radiometer is modeled by the appropriate radiative transfer equation. For an airborne, downward looking radiometer, the measured emission in a scattering-free medium will be due to the surface emission, which is attenuated by the atmosphere, the upward atmospheric emission, and the downward atmospheric and cosmic background emission that is reflected by the surface in the viewing direction of the radiometer’s antenna. The microwave radiative transfer equation has the form,

 
formula

where Tsfc is the physical temperature of the surface (K), TC is the microwave cosmic background temperature (∼2.7 K), ɛ is the surface emissivity, and θ is the incidence angle of the radiometer (Ulaby et al. 1981). Specular surface reflection is assumed in (4). The upwelling and downwelling (TBUP and TBDOWN) atmospheric brightness temperature can be calculated given the temperature and extinction/absorption coefficient profile. The optical depth, τ, is the integral of the total extinction coefficient along a path, where the total extinction coefficient is the sum of the contributions from precipitation, cloud liquid water and atmospheric gasses. The absorption coefficient is related to the extinction coefficient by κabs = (1 − ω0)κext, where ω0 is the single scattering albedo. Equation (4) neglects the effects of scattering, with the exception of the scattering coefficient in the extinction term. The single scattering albedo at 10.7 GHz is less than 0.08 for rain rates less than 25 mm h−1, which is why scattering can be neglected with little error. A possible situation in which this assumption would not be valid is the case of a large volume of ice particles that scatter primarily in the forward direction. The brightness temperature in this case would be underestimated by Eq. (4).

The drop size distribution at each range gate determines the volume extinction and absorption coefficient at 10.7 GHz due to the rain according to

 
formula

assuming the proper Mie absorption and extinction efficiency based on the thermodynamic makeup of the particles. Profiles of the extinction and absorption coefficient from the precipitation, derived from the PR-2 DSD retrieval, are summed with the volume extinction coefficient from the assumed cloud liquid water and gaseous extinction coefficient to determine the upwelling and downwelling atmospheric brightness temperatures and the optical depth. The radiative transfer equation, (4), is inverted to estimate the surface emissivity, which has monotonic dependence on the surface wind speed. In practice, the wind speed is directly estimated using the iterative inversion method described in section 3.

The wind-roughened sea surface emissivity model used in this study is from Pandey and Kakar (1982). For a horizontally polarized wave, the increase above the specular case in the sea surface emissivity will be due to the surface roughness from the waves and patches of sea foam. The empirical surface emissivity model is described by

 
formula

where Δeroughness is the additional surface emissivity introduced by wind driven waves, F is the fractional sea foam coverage, and efoam is the emissivity of sea foam. Sea foam has an emissivity near unity in the microwave. Therefore, the extent to which the footprint is covered by foam will have a large impact on the emissivity of the observed area. Sea foam is generally present when the wind speeds become larger than 7 m s−1 and the fractional foam coverage increases monotonically with wind speed due to breaking waves and roughening of the surface. In Pandey’s model, the expression for the surface roughness term is described by an empirical expression from Hollinger (1970) and the expression for the emissivity of the sea foam is from Stogryn (1972). Seasat–SMMR (Scanning Multichannel Microwave Radiometer) brightness temperature (TB) measurements were used to determine a quadratic fit with wind speed for the foam fraction.

Because LRR oversamples each pixel and has a higher resolution than PR-2, there are many LRR points at slightly different incidence angles for each atmospheric attenuation profile derived from the PR-2. The inversion seeks to find the near surface wind speed that minimizes the RMS difference between the modeled and measured TB for all data points falling within the volume corresponding to the attenuation profile. Because the atmospheric optical depth and upwelling and downwelling brightness temperatures are determined for the slant path corresponding to the PR-2 incidence angle, these values are scaled by secθ for the slightly different incidence angles corresponding to the LRR measurements. The LRR incidence angles are on average ±1° from the PR-2 incidence angle. Convergence is reached when the RMS difference between the modeled and measured TBs is less than 0.5 K. The maximum number of iterations allowed is 13. For cases where the atmospheric absorption is overestimated and the value of the wind speed goes negative, the retrieval is set to 0 m s−1.

A possible source of error arises from the plane parallel assumption when retrieving the wind speed under isolated convective cells at high incidence angles, toward the edge of the scan. In this case, especially near the leading or trailing edge of the storm in the direction orthogonal to the DC-8 flight line, there could be a large difference in the magnitude of the upwelling and downwelling brightness temperatures. A simple first-order correction can be used, in which the downwelling brightness temperature is set to that downwelling brightness temperature determined from PR-2 at the adjacent pixel in the scan with the higher incidence angle. Errors arising from the plane parallel assumption should be minimal for incidence angles close to nadir and will be minimal for most incidence angles when the spatial inhomogeneity of the rain is low, such as in stratiform rain.

5. Melting layer model

Determining the absorption and scattering parameters in the melting layer presents a more complicated situation as compared to the snow or rain layer. The melting layer generally forms in stratiform situations when the vertical and horizontal air velocities are small. This allows a layer of melting particles to exist just below the height of the freezing level. The melting layer is generally 300–600 m deep. The effect of the melting layer on radar backscattering was first observed as a relatively bright layer on the radar operators analog display tube, giving it the name bright band (Rinehart 1997). This bright layer is caused by increased backscatter from the melting particles. As a low-density, large-diameter snow particle begins to melt, it takes on electromagnetic properties closer to that of water, while still maintaining a large diameter relative to a rain drop of the same mass. In this way, the melting particles act like abnormally large rain drops, causing the melting layer to have increased scattering and absorption. To determine the absorption and scattering in the melting layer, one must use a thermodynamic model to determine the melted mass fraction as a function of height below the freezing level and a dielectric mixing formula to determine the absorption and scattering coefficients for a given melting particle. The thermodynamic formula will generally drive the thickness of the melting layer and the choice of the dielectric model will generally determine the intensity of the absorption and scattering. The components of the melting layer model have been well described in previous work and will only be summarized here (Bauer et al. 1999; Fabry and Szyrmer 1999; Szyrmer and Zawadzki 1999; Olson et al. 2001; Battaglia et al. 2003).

a. Melting layer thermodynamic model

The melting rate of a snow particle that falls below the 0°C isotherm can be described in terms of a thermodynamic heating budget (Mitra et al. 1990). Melting occurs so long as the sensible heat transfer at the particle surface is greater than the evaporative cooling. This model provides the melted mass fraction of a given diameter melting snowflake as a function of height below the 0°C isotherm (Willis and Heymsfield 1989; Mitra et al. 1990; Battaglia et al. 2003). The melting rate for a given diameter snowflake is determined by its initial snow density, fall velocity, the ventilation coefficient and the environmental conditions, such as the temperature lapse rate and relative humidity.

The ventilation coefficient describes the enhanced vapor and sensible heat gradients experienced by a particle moving relative to air as compared to a still one (Wusheng and Wang, 1999). Several parameterizations have been used by various authors and currently there is no experimental evidence to favor one over another. Three such parameterizations that are tested in this study are from Rutledge and Hobbs (1983), Mitra et al. (1990), and Szyrmer and Zawadzki (1999).

Additionally, there exist many parameterizations for snowflake density based on in situ observations, which is reflective of the numerous ice habits that a frozen particle can take. Typically these models are of the form

 
formula

where Ds is the snowflake diameter. Model parameters from several investigators are shown in Table 1. A lower density snowflake will melt faster than a higher density snowflake of the same diameter. Therefore, dense snowflakes will produce a deeper melting layer as compared to fluffy aggregates. For a given particle diameter the Magono and Nakamura (1965) model produces higher density snow flakes and the University of Wisconsin–Nonhydrostatic Modeling System (UW–NMS) model produces the least dense snow flakes. The other models are in between. The value for the snow density is truncated at 0.9 and 0.005 g cm−3 as suggested by Klaassen (1988).

Table 1.

Snow density parameters for Eq. (7). The units for snow flake diameter are in cm.

Snow density parameters for Eq. (7). The units for snow flake diameter are in cm.
Snow density parameters for Eq. (7). The units for snow flake diameter are in cm.

Battaglia et al. (2003) provides a parameterization for the velocity of a melting particle derived from Mitra et al.’s (1990) experimental work and shows that the output of the melting layer model is not very sensitive to the parameterization for the fall velocity. The melting rate is also affected by the environmental parameters. A high lapse rate will cause the particle to melt faster. Conversely, if the relative humidity is low, the particle melting rate will slow because of the sensible heat flux being balanced by the evaporative cooling from the meltwater. The fall velocity parameterization and the environmental parameters, which are determined from a nearby radiosonde sounding, are fixed in the following analysis.

b. Melting layer electromagnetic models

Determining the absorption and scattering coefficients for a melting particle requires that one define an effective dielectric constant for the mixture of air, ice, and water. The thermodynamic formula will give the melting particle diameter and volume fraction of the three constituents at any stage in the melting process and a proper dielectric mixing formula must be applied to determine the dielectric constant of the particle. There are several choices of dielectric mixing formulas in the literature. The most commonly used form to model melting particles is the Maxwell–Garnet (1904, hereafter MG) formula, which was generalized for elliptical inclusions by Bohren and Battan (1982). In the MG formula, the choice of which component is the matrix and which is the inclusion will affect the resulting dielectric of the mixture. For a three-component mixture, such as a melting particle, the MG formula must be applied twice. In all, there are 12 different ways to formulate the MG equation to accommodate air, ice, and water. Some commonly used combinations are summarized in Table 2 and given common notation that will remain consistent throughout the rest of this study (Meneghini et al. 1997; Bauer et al. 1999; Olson et al. 2001; Battaglia et al. 2003). The combinations in which water forms the matrix for the final application of the formula will produce the stronger absorption and scattering properties compared to the combinations in which ice or air make up the matrix.

Table 2.

Maxwell–Garnet combinations used in this study.

Maxwell–Garnet combinations used in this study.
Maxwell–Garnet combinations used in this study.

Two other models that appear in the literature are the Meneghini and Liao models (Meneghini and Liao 1996, hereafter ML96) and (Meneghini and Liao 2000, hereafter ML00), and the Fabry and Szyrmer (1999) coated sphere model (FS). The Meneghini and Liao models are not used in this study because they have been shown to underestimate the intensity of both the reflectivity peak and total attenuation in the melting layer (Olson et al. 2001; Battaglia et al. 2003). The Fabry and Szyrmer model is a reasonable candidate, as it has compared favorably with measured reflectivity profiles in the melting layer from a 0.915-GHz radar profiler (Battaglia et al. 2003) and radar derived estimates of the Path Integrated Attenuation (PIA) at 13.8 GHz (Olson et al. 2001). The model formulation is summarized below.

The density of a snowflake is generally assumed to decrease with particle diameter. To model a snow particle whose density varies throughout its diameter, Fabry and Szyrmer (1999) divided the particle into a core and a shell with different densities, where the density transition takes place at αsDs for a snow particle of diameter Ds. Through mass conservation, the density of the snowflake inner core ρins and the outer shell ρouts can be determined from the total snow particle density ρs. The initial location of the density transition, αs, was arbitrarily set to 0.5 in this study. During the melting process, the density of the inner core and the outer shell will change as the melted mass fraction of the particle changes. In this model, it is assumed that the meltwater fraction, fw, is the same in the core and the shell. The density of the inner and outer core of the melting particle can be expressed as a function of the initial snowflake core and shell densities and the melt water fraction assuming mass conservation. This gives a unique volume fraction of air, ice, and water for the core and the shell, which can be used to determine the effective dielectric of each. The dielectric mixing formula for a melting particle consisting of a core with dielectric ɛcore and a shell with dielectric ɛshell is then used to determine the effective dielectric of the whole particle. The dielectric of the core and the shell can be determined using one of the MG combinations (e.g., MG1, MG2, MG3, etc.), or any other suitable mixing formula. This model provides a means to blend two dielectric models that individually over and underestimate the scattering, to provide, in a reasonable way, a better fit to radar data in the melting layer.

c. Fitting melting layer models to PR-2 profiles

The melting layer typically contributes a majority of the total atmospheric attenuation in light stratiform rain where the freezing level is close to the surface. The melting layer models discussed previously can produce vastly different amounts of absorption and scattering. Therefore, to produce reasonable wind speed and rain rate estimates, it is necessary to find a melting layer model that predicts a realistic extinction coefficient profile. Because there are no direct measurements of the extinction coefficient itself, the model chosen is that which best fits the measured backscattering profiles at 13.4 and 35.6 GHz. To constrain the relationship between extinction and backscatter, the DSD must be known. The DSD in the melting layer can be determined by assuming stationary conditions; that is, constant inflow and outflow of mass at each level in the melting layer. This is equivalent to the condition

 
formula

Thus, given the parameterized drop size distribution at the base of the melting layer, the DSD at any height in the melting layer is determined through (8) using the expressions given for the velocity of melting and melted particles. The thermodynamic model is used to determine the volume fraction of air, ice and water that a particle will have as a function of its diameter and distance below the freezing level. These values are calculated for liquid equivalent drop diameters of 0.02 to 10 mm in 0.02-mm increments and in 10-m steps of height below the freezing level from 0 to 1000 m. A dielectric mixing model is then used to determine the Mie efficiencies for the melting particles at 10.7, 13.4, and 35.6 GHz.

To fit the models to the PR-2 measured profiles, the DSD can be determined at the base of the melting layer through the dual-frequency method described in section 3, but the reflectivities must first be corrected for attenuation by the melting and snow layers. This is accomplished by first estimating D0 using only the measured reflectivity at 13.4 GHz with an average value of N*0 = 5.0 × 106 m−4. The same iterative inversion method is used in the single frequency DSD retrieval as in the dual-frequency DSD retrieval with the exception that only one parameter is estimated. The initial DSD is then used with the particular melting layer model to obtain an estimate of the total attenuation in the melting layer at both PR-2 frequencies. These values are used to correct the basal level reflectivities, which are then used to estimate N*0 and D0 with the dual-frequency version of the inversion algorithm. This is the final DSD that is used as the input to the melting layer model. It was observed that iterating this process again does not appreciably change the estimates of N*0 and D0 at the base of the melting layer, suggesting that the fitting procedure has converged. An example of this fitting procedure is shown in Fig. 1 for both frequencies. In Fig. 1, measured PR-2 profiles are shown with the reflectivity profiles predicted using the MG formulas listed in Table 2 with the Magono and Nakamura snow density model and the Mitra ventilation coefficient. The temperature lapse rate was assumed to be 7.7 K km−1, which was determined from the coincident radiosonde observation (RAOB), and the relative humidity was assumed to be 100%. Figure 2 shows an example of the corresponding extinction coefficient profile in the melting layer predicted by the models in Fig. 1 at 10.7 GHz. A strong increase of the extinction coefficient in the melting layer is observed at all frequencies. The freezing level is estimated to be at 1503 m for this profile. It is observed that the dielectric formulas in which water is included in the matrix for the final application of the formula (MG3 and MG4) produce the strongest extinction peaks. The dielectric models that produce the strongest reflectivity peaks also produce the strongest extinction peaks.

Fig. 1.

Example of the melting layer model fits to the PR-2 data at (left) 13.4 and (right) 35.6 GHz. Shown are the model fits using the MG models MG1, MG2, MG3, and MG4 described in Table 2. The PR-2 data are the lines with the dots.

Fig. 1.

Example of the melting layer model fits to the PR-2 data at (left) 13.4 and (right) 35.6 GHz. Shown are the model fits using the MG models MG1, MG2, MG3, and MG4 described in Table 2. The PR-2 data are the lines with the dots.

Fig. 2.

The 10.7-GHz extinction coefficient profile in the melting layer for the same PR-2 profile and MG models in Fig. 1.

Fig. 2.

The 10.7-GHz extinction coefficient profile in the melting layer for the same PR-2 profile and MG models in Fig. 1.

To assess how well each melting layer model fits the PR-2 data, two key features are analyzed. First, the melting layer model should be able to reproduce the peak of the scattering at both frequencies. This is mainly driven by the choice of the dielectric model. How well each model fits the peak will be assessed by comparing the maximum reflectivity from the model, regardless of where it occurs, to the maximum reflectivity from the measured profile. The second feature that is assessed is the width of the model relative to the measurements. This will be mainly driven by the choice of the ventilation coefficient and the initial snow density, as well as the environmental parameters (i.e., lapse rate and relative humidity). To assess how well each model fits the measured width, the height difference at the half maximum point is compared. Here the half maximum refers to half of the difference between the basal reflectivity and the reflectivity peak. Both of these features are important in order to model the attenuation in the melting layer. For example, a particular model could exactly fit the scattering peak, but have a much larger width compared to the measured profile. This would produce an overestimate of the total attenuation in the melting layer, even though the peak of the extinction was correct. These two metrics were chosen instead of the RMS difference between the modeled and measured profiles because the RMS will be sensitive to errors in the determination of the melting layer boundaries.

d. Results of melting model comparisons

The largest uncertainties in the melting layer model are the parameterization for the ventilation coefficient, the initial snow density and the dielectric mixing model to describe the absorption and scattering properties of the partially melted hydrometeors. The choice of the dielectric formula drives the strength of the absorption and scattering for the melting layer model, as illustrated in Figs. 1 and 2. The choice of the ventilation coefficient and the initial snow density are second-order terms.

An initial analysis is conducted with a broad field of candidates—31 different combinations representing five dielectric mixing formulas. The analysis is conducted on approximately 100 adjacent stratiform profiles with basal reflectivities varying over 25–31 dBZ. The lapse rate was assumed to be 7.7 K km−1 and the air was assumed to be saturated. Table 3 provides a summary of the results for the different dielectric models using the Mitra ventilation coefficient and the Magono and Nakamura snow density model. In the table, the peak bias is the model reflectivity peak minus the PR-2 reflectivity peak. The fraction of the opacity is the fraction of the total atmospheric extinction at 10.7 GHz due to the melting layer. The mean wind speed is the average retrieved surface wind speed from the 100 profiles. It can be seen that the MG4 and MG3 models overestimate the reflectivity peak at both frequencies, by 4 and 2.6 dB at 13.4 GHz, respectively. These models are also observed to overestimate the total atmospheric absorption at 10.7 GHz, which is evident in the nonrealistic wind speed retrievals (negative wind speeds were set to zero). The MG2 and MG1 models underestimate the reflectivity peak at both frequencies. These models also produce the least amount of attenuation in the melting layer. Of the dielectric models shown in Table 3, the FS core–shell model produces the best agreement with the observed reflectivity peaks at both frequencies. It is biased low by 1.3 dB at 13.4 GHz and 0.4 dB at 35.6 GHz. The current configuration has a strongly reflecting core, MG3, surrounded by a weakly reflecting shell, MG2. The reflectivity peak predicted using the FS model is between the reflectivity peaks predicted individually by the MG dielectric formulas of the core and the shell. From this initial analysis, it can be inferred that a method that blends two Maxwell–Garnett formulas will produce results that best fit the measured profiles, of which the Fabry–Szyrmer core–shell model is an example. Changing the location of the density transition in the snow particle, which was arbitrarily set by Fabry–Szyrmer to 0.5, will modify the relative contributions from the core and the shell. Because the current configuration has a strongly reflecting core, if the initial snow density transition is increased, the scattering and absorption properties of the particle will increase. It was determined that changing αS to 0.75 provided a better fit to the PR-2 data in terms of the metrics discussed above. This model will henceforth be referred to as the FS-a75 model. It is not unreasonable to expect the density transition of the melting particle to occur toward the edge of the particle. Observations of melting snow flakes in a wind tunnel suggest that the melt water first forms on the tips of the ice crystals at the periphery of the flake (Mitra et al. 1990). As melting continues, the melt water is collected, through capillary forces, at the center of the flake around the ice crystal structure. Melting proceeds in this fashion until the crystal supports collapse and the particle forms a raindrop.

Table 3.

Fit to PR-2 data for different dielectric models.

Fit to PR-2 data for different dielectric models.
Fit to PR-2 data for different dielectric models.

A secondary analysis was conducted with the FS-a75 model, using three parameterizations for the ventilation coefficient and three initial snow density models, on 525 adjacent stratiform reflectivity profiles. Table 4 shows the results using the Barthazy initial snow density for the Mitra, Rutledge, and Szyrmer and Zawadzki (SZ) ventilation coefficients. Likewise, Table 5 shows the results using the Locatelli initial snow density and Table 6 shows the results using the Magono initial snow density. Of the three density models, the Magono model will produce overall the densest snow flakes, followed by the Locatelli model and then the Barthazy model. The SZ ventilation coefficient produces a melting layer, whose depth is almost independent of snow density, which is evident in the tables as the width bias being nearly the same for the three density models. The low-density Barthazy and Locatelli models tend to underestimate the width of the melting layer by about 70 m using the Mitra ventilation coefficient and 100 m using the Rutledge ventilation coefficient. The Magono model produces a deeper melting layer due to the slower melting rate of the denser particles. The width using the Magono density model and the Mitra ventilation coefficient is slightly underestimated by 48 m. The Mitra and the Rutledge ventilation coefficients produce similar results for the reflectivity peak at both frequencies for a given snow density model, which is evident in the first two columns of each table. However, they do not produce similar melting layer depths. The Rutledge ventilation coefficient consistently produces a shallower melting layer as compared to the Mitra ventilation coefficient. The Barthazy snow density model overestimates the 13.4-GHz brightband peak by 0.73 dB using the Mitra ventilation coefficient and 0.78 dB using the Rutledge ventilation coefficient and it underestimates the 35.6-GHz brightband peak by 0.91 dB using the Mitra ventilation coefficient and 0.83 dB using the Rutledge ventilation coefficient. The Locatelli snow density model is better able to fit the 35.6 GHz scattering peak, but overestimates the 13.4-GHz scattering peak slightly more than with the Barthazy density model. The SZ ventilation coefficient with the Barthazy and Locatelli snow density models is biased by less than 0.1 dB from the measured 13.4-GHz reflectivity peak, but it underestimates the 35.6-GHz reflectivity peak for both snow density models. Notice in Tables 3 and 4 that even though the predicted 13.4-GHz reflectivity peak using the SZ parameterization is almost 1 dB lower than that using the Mitra and Rutledge parameterizations, the total attenuation in the melting layer is higher using the SZ parameterization because of the larger depth of the melting layer compared to the other two. The Magono snow density model performs the best overall. The peak bias at 13.4 and 35.6 GHz is low for the Mitra and the Rutledge parameterizations. The 13.4-GHz bias is −0.16 dB with the Mitra ventilation coefficient and −0.13 dB using the Rutledge ventilation coefficient and the 35.6-GHz bias is 0.33 dB using the Mitra ventilation coefficient and 0.41 dB using the Rutledge parameterization. The reflectivity peak is underestimated by 1.2 dB at 13.4 GHz and 0.34 dB at 35.6 GHz for the SZ parameterization. The combination that performs the best, based on fitting the width and the peak of the measured reflectivity profiles, is the FS-a75 model with the Magono snow density model and the Mitra ventilation coefficient. It should be noted that the width of the melting layer is also dependent on the assumed lapse rate and relative humidity profile. The peak predicted by the model has little sensitivity to the environmental variables. Figure 3 shows a plot of the brightband peak versus the basal reflectivity for this model and for the observations. The RMS difference of the reflectivity peaks between the model and the measurements is 1.6 dB at 13.4 GHz and 1.1 dB at 35.6 GHz. For the most part, the model is able to predict the reflectivity peak at both frequencies over a range of basal reflectivities. The scatter of the observations is due to measurement noise and departure from the stationary and mass conservation assumptions. Notice that there is a group of observations around ZBA = 20 dBZ that do not fit the general trend. That is, the difference between the scattering peak and the basal reflectivity is much higher than for other observations at this basal level and higher than predicted by the model. This is most likely because the stationary assumption is not valid at this location due to a perturbation from the wind field. Overall, the reflectivity peak predicted by the model is in agreement with the observations from PR-2 at both frequencies.

Table 4.

Metrics for the FS-a75 model with the Barthazy initial snow density.

Metrics for the FS-a75 model with the Barthazy initial snow density.
Metrics for the FS-a75 model with the Barthazy initial snow density.
Table 5.

Metrics for the FS-a75 model with the Locatelli initial snow density.

Metrics for the FS-a75 model with the Locatelli initial snow density.
Metrics for the FS-a75 model with the Locatelli initial snow density.
Table 6.

Metrics for the FS-a75 model with the Magono initial snow density.

Metrics for the FS-a75 model with the Magono initial snow density.
Metrics for the FS-a75 model with the Magono initial snow density.
Fig. 3.

Modeled (*) and observed (x) brightband peak reflectivity versus observed basal reflectivity for (left) 13.4 GHz and (right) 35.6 GHz. The mean and RMS of the difference is −0.16 and 1.6 dBZ at 13.4 GHz and 0.33 and 1.1 dBZ at 35.6 GHz. The slope and intercept for the model are 1.0 and 9.6 at 13.4 GHz and 0.70 and 9.3 at 35.6 GHz.

Fig. 3.

Modeled (*) and observed (x) brightband peak reflectivity versus observed basal reflectivity for (left) 13.4 GHz and (right) 35.6 GHz. The mean and RMS of the difference is −0.16 and 1.6 dBZ at 13.4 GHz and 0.33 and 1.1 dBZ at 35.6 GHz. The slope and intercept for the model are 1.0 and 9.6 at 13.4 GHz and 0.70 and 9.3 at 35.6 GHz.

6. Case study

A case study is presented to demonstrate the retrieval algorithm. The case study is for a region of stratiform rain with an embedded convective cell. A plan view of the reflectivities around 700 m and a vertical cross section at y-index 13 (referenced to the y axis in the plan view) are shown in Fig. 4. The storm top (height of the first detectable return) is at approximately 4 km above the surface, although the optical cloud top may be well above this level. The horizontal extent of the storm is approximately 60 km, from index 50 to index 200 in Fig. 4. Each index is roughly 450 m. There is an embedded convective cell with higher reflectivities around index 140. There are rain-free regions on either side of the storm, although these areas most likely contain variable amounts of cloud liquid water. The rain-free regions can be used as a metric to assess the performance of the wind speed retrievals under the rain. That is, we should not expect a large jump in the retrieved wind speed between the raining and nonraining portions of the image. Some variability in the retrieved wind speed is expected due to variations in the amount of cloud liquid water (CLW), which is assumed to be constant across the image. For these retrievals, the gaseous absorption profile is determined from a coincident RAOB sounding. A liquid cloud with a constant water density of 0.1 g m−3 is assumed from 500 to 1500 m and an ice cloud with the same density from 1500 to 6000 m.

Fig. 4.

(top) Plan view and (bottom) vertical cross section of ray 13 for the PR-2 13.4-GHz reflectivities. Each index is roughly 450 m.

Fig. 4.

(top) Plan view and (bottom) vertical cross section of ray 13 for the PR-2 13.4-GHz reflectivities. Each index is roughly 450 m.

Figure 5 shows a plot of the average retrieved rain rate in the rain layer for y-index 13 (nadir observation). The rain rates in the stratiform region are less than 1 mm h−1. The rain rate peaks at about 12 mm h−1 in the embedded convective cell. Figure 6 shows a plot of the wind speed retrieval under the rain. The dashed line is the wind speed retrieval assuming a nonprecipitating atmosphere; that is, assuming only absorption from oxygen and water vapor and no correction for the presence of rain. The thick line is the wind speed retrieval that has been corrected for the presence of the rain, snow, and melting layers. Overall, the retrieval algorithm produces wind speed retrievals under the rain that are physically reasonable, as compared to the rain-free regions on either side of the storm. The retrieved wind speed is fairly constant in the transition from the stratiform region to the embedded convective cell, even though the rain rate is changing from 1 to 12 mm h−1. The wind retrieval peaks near index 75 and 175 could be due to residual atmospheric brightness that is not accounted for. Figure 7 shows the average retrieved median drop diameter in the rain layer and Fig. 8 shows the average retrieved N*0 in the rain layer. The transition from the stratiform rain to the convective rain is marked by an increase in N*0 and in D0. The increase in N*0 is forced due to the different a priori bounds on N*0 for convective rain. The median drop diameter is about 0.09 cm in the stratiform rain and increases to about 0.12 cm in the convective cell. N*0 is around 4.46 m−4 in the stratiform region and about 18.56 m−4 in the convective region. It was observed that the RMS error in the DSD retrieval at all heights in the convective rain layer increased when the stratiform a priori bounds (lower N*0) were used, justifying the increase in N*0.

Fig. 5.

Rain layer averaged retrieved rain rate for PR-2 ray 13.

Fig. 5.

Rain layer averaged retrieved rain rate for PR-2 ray 13.

Fig. 6.

The retrieved wind speed for ray 13 is shown as the black line. The dotted line shows the wind speed retrieval without correction for precipitation.

Fig. 6.

The retrieved wind speed for ray 13 is shown as the black line. The dotted line shows the wind speed retrieval without correction for precipitation.

Fig. 7.

Rain layer averaged D0 retrieved for ray 13.

Fig. 7.

Rain layer averaged D0 retrieved for ray 13.

Fig. 8.

Rain layer averaged N*0 retrieved for ray 13.

Fig. 8.

Rain layer averaged N*0 retrieved for ray 13.

To illustrate the effect of the melting layer on the wind speed retrieval, a retrieval was performed with the melting layer removed. The precipitation profile is assumed to transition directly from rain to snow at the base of the melting layer. Figure 9 shows the difference between the wind speed retrieval with and without the melting layer present. The difference is also shown when precipitation is completely removed from the algorithm (i.e., assume clear-sky conditions). Notice that the wind speed retrieval is driven by the absorption in the melting layer. Only about 1 m s−1 of difference is explained by the rain (the difference between the two lines), whereas 2–4 m s−1 of difference is explained by the melting layer.

Fig. 9.

Difference in wind speed retrieval when the melting layer is removed and when both the melting layer and precipitation are removed.

Fig. 9.

Difference in wind speed retrieval when the melting layer is removed and when both the melting layer and precipitation are removed.

A cross -track image of the retrieved rain rate, averaged in the rain layer, is shown in Fig. 10. The rain rates are very light in the stratiform region, generally less than 1 mm h−1 across the swath. The rain rates peak around 15 mm h−1 in the embedded convective cell. The retrieved wind speed image under the rain is shown in Fig. 11. There are a few pixels in the convective region where an anomalously low value of the surface wind speed was retrieved, although these represent only a small portion of the retrieved pixels. The surface wind speeds are observed to increase radially outward around the convective cell on the northeast side of the storm, from x, y indices (160, 1) to (138, 23). This could be due to a gust front propagating outward from the convective core. The propagation of the front would be strongest on the northeast side because this is the direction of the environmental wind field.

Fig. 10.

Cross-track image of the retrieved rain rate averaged in the rain layer.

Fig. 10.

Cross-track image of the retrieved rain rate averaged in the rain layer.

Fig. 11.

Cross-track image of the retrieved surface wind speed.

Fig. 11.

Cross-track image of the retrieved surface wind speed.

A scatterplot of the retrieved wind speeds with and without the precipitation correction is shown in Fig. 12 as a function of the peak reflectivity above 900 m. For the stratiform case, this represents the brightband peak. The uncorrected wind speeds are highly correlated with the reflectivity, and rise to about 28 m s−1. The retrieved wind speeds that have been corrected for precipitation are uncorrelated with the reflectivity and are steady around 9 m s−1. The mean and standard deviation of the retrieved wind speeds that have been corrected for precipitation are shown in Fig. 13 in 1-dBZ bins of peak reflectivity from 20 to 46 dBZ. The standard deviations are around 2 m s−1 below 35 dBZ, and around 3 m s−1 above. Some of this variance is due to uncertainties introduced by the retrieval algorithm and some is due to the natural variation of the winds. The fact that retrieved wind speed is fairly constant even though the peak reflectivity is changing from 20 to 46 dBZ provides confidence that the retrieval algorithm is properly accounting for the attenuation from precipitation. This implies that the melting layer model, which was chosen based on its fit to the PR-2 reflectivities, provides a reasonable estimate for the attenuation at 10.7 GHz. Also, this provides confidence that the rain retrievals in the convective area are reasonable, since there is no noticeable bias when the reflectivities are greater than 40 dBZ.

Fig. 12.

Retrieved wind speeds with (dot) and without (x) correction for precipitation.

Fig. 12.

Retrieved wind speeds with (dot) and without (x) correction for precipitation.

Fig. 13.

Mean and std dev of the retrieved wind speed grouped into 1-dBZ bins of peak reflectivity.

Fig. 13.

Mean and std dev of the retrieved wind speed grouped into 1-dBZ bins of peak reflectivity.

7. Sensitivity of retrievals

The sensitivity of the retrievals to measurement noise, the assumed cloud liquid water, sea surface temperature, and freezing-level height was assessed using a Monte Carlo technique. The retrieved average rain rate varied by about 6% when random additive Gaussian noise with a standard deviation of 0.5 dB was added to the reflectivities. The uncertainty introduced into the wind speed retrievals from the radar noise was determined to be 0.7 m s−1 for light stratiform rain and 1.7 m s−1 for moderate convective rain. The uncertainty in the wind speed in light stratiform rain has an almost perfect negative correlation with the uncertainty in the attenuation from the melting layer and is almost uncorrelated with the uncertainty in the rain retrieval, suggesting that radiometric retrievals in light stratiform precipitation are almost completely driven by the model for the melting layer. This is discussed in more detail in section 8, as a technique to constrain the melting layer if the wind speed is known.

It was found that the sensitivity of the wind speed retrievals to uncertainty in the assumed SST and the radiometer NEΔT was minimal. The uncertainty of the freezing-level height determined by the brightband detection algorithm was estimated to be 50 m in stratiform rain and 100 m in convective rain. This uncertainty introduced a 0.4% error in the average rain retrieval for stratiform rain and a 2.1% error in the average rain retrieval for convective rain. Uncertainties in the freezing level introduced a 0.75 m s−1 uncertainty in the wind speed retrieval (at 7 m s−1) in stratiform rain and a 2.4 m s−1 uncertainty in the wind speed retrieval in moderate convective rain.

The retrievals could contain errors that are spatially correlated because of the uncertainty in the cloud liquid water content in the profile. If the environmental CLW is higher than the value assumed in the retrieval algorithm, then the rain estimates will be consistently biased low and the wind speed estimates will be consistently biased high. The wind speed bias for a given amount of CLW will be larger for higher rain rates. The potential error sources for the rain rate and wind speed retrievals are summarized in Tables 7 and 8. The largest source of random error in the rain rate is due to the noise in the reflectivities. The rain rates will be biased high if the assumed CLW is biased high and vice versa if the CLW is biased low.

Table 7.

Estimated uncertainty in average rain retrieval for different error sources for the top-down approach.

Estimated uncertainty in average rain retrieval for different error sources for the top-down approach.
Estimated uncertainty in average rain retrieval for different error sources for the top-down approach.
Table 8.

Uncertainty in the wind speed retrieval due to different error sources for the top-down approach.

Uncertainty in the wind speed retrieval due to different error sources for the top-down approach.
Uncertainty in the wind speed retrieval due to different error sources for the top-down approach.

The uncertainty in the wind speed retrieval due to the error sources investigated is summarized in Table 8 at a wind speed of 7 m s−1 for rain rates of 0.84 and 14.3 mm h−1. The largest source of error is the uncertainty in the CLW. Fortunately, this uncertainty can be minimized by the addition of another radiometer channel. The second largest source of error for the convective profile is the uncertainty in the height of the freezing level. For the stratiform profile, the uncertainty in the height of the freezing level is almost the same as that caused by reflectivity noise. Uncertainty in the SST and the radiometer NEΔT are minimal sources of error.

The assumptions made concerning the DSD will affect the rain rate and wind speed retrievals. Dual-frequency reflectivity measurements are used to determine two of the three free parameters of the normalized gamma DSD. This requires an assumption for the value of the third parameter. In this study, the dispersion factor is set to a value of 0.99 that is constant with height and location. To assess the sensitivity of the retrievals to differing values of μ, 75 realizations of the retrieval were performed on nine convective profiles with a normally distributed random value for μ with a mean of 0.99 and standard deviation of 1.0. The resulting average standard deviation of the retrieved rain rate was 13% and the average standard deviation of the retrieved wind speed at 7 m s−1 was 2.7 m s−1. This suggests that large variations in μ can cause large uncertainties in the retrievals. Quantifying the amounts by which μ varies for differing storm conditions is important to assess the impact of this uncertainty and if a third independent measurement of the precipitation volume is required.

To assess how much of the observed variation of the wind speed retrievals can be explained by the aforementioned errors, a simple error budget is constructed. The standard deviation of the wind speed retrievals from the case study was approximately 2 m s−1 in the stratiform rain (see Fig. 13). One would expect little variation of the true wind speed in the stratiform region, suggesting that most of the 2 m s−1 of standard deviation is due to error. The root-sum-of-squares (RSS) of the estimated random uncertainties in the wind speed for stratiform rain (see Table 8); that is, the instrument noise and the freezing level height, is 1.1 m s−1. A small variation in the environmental CLW, ±0.2 g m−3 in addition to the random uncertainty, would explain the 2 m s−1 of standard deviation observed, assuming no variation in the actual wind field. This value for CLW fluctuation is not unreasonable, especially at the small spatial scale of observation. An image of retrieved CLW from the Advanced Microwave Scanning Radiometer (AMSR-E) showed pixel-by-pixel variations of integrated CLW from 0 to 0.5 mm at a spatial resolution of 30 km in the vicinity of the retrieval area at nearly the same time (corresponding to 0–0.5 g m−3 for a uniform 1 km thick liquid cloud). The variations of CLW observable at a spatial resolution of 400 m are most likely greater than this.

The standard deviation of the retrieved wind speed from the case study in the convective rain was 3 to 4 m s−1. Near convective rain, one would expect some variation in the surface wind speed, so some of the deviation from the mean may be due to actual variations in the wind field, making the construction of an error budget more ambiguous. The RSS of the random uncertainties for the convective values shown in Table 8 is 2.9 m s−1. The remainder of the observed 3–4 m s−1 of variation is most likely explained by actual variations in the wind field and variations in the CLW.

8. Using the algorithm to constrain the melting layer model

In addition to the variables discussed in the previous section, the retrieved wind speed is also very sensitive to the choice of the melting layer model (see details in section 5). Actually, the uncertainty on the wind speed retrieval from the melting layer model trumps all other error sources explored in section 7 for light to moderate stratiform rain. This presents a robust technique for constraining the melting layer using the retrieval algorithm if the surface wind field is known to a certain accuracy. The melting layer model comparisons from section 5 are used to estimate the sensitivity of the derived wind speed to the 10.7-GHz attenuation in the melting layer (note the nearly perfect negative correlation between melting layer opacity and wind speed in Tables 3 –6). The average sensitivity value for the 525 profiles tested is determined to be 0.0024 Np m s−1 for wind speeds near 8 m s−1 and light to moderate rain rates. The average 10.7 GHz melting layer optical depth of these profiles for the selected FS-a75 model is 0.0112 Np. The average optical depths predicted using the other ventilation coefficient and snow density combinations with the FS-a75 model range from 0.0101 to 0.0175. This means that if the surface wind speed is known to 1 m s−1 (from a coincident buoy or other means), only those models that predict attenuation values within about 20% of each other would be ambiguous. This is illustrated graphically in Fig. 14. Here the melting layer optical depth is averaged in 1-dBZ bins of basal reflectivity and plotted for each combination listed in Tables 4 –6. The sensitivity of the optical depth to wind speed is determined for each basal reflectivity bin. The selected model is shown as a bold line and the errors bars placed on this line represent the ambiguous region for a wind speed knowledge of 1 m s−1. Using this example, only four combinations would be ambiguous with the wind speed constraint. This constraint, coupled with the dual-frequency reflectivity profile constraints, present a powerful method for distinguishing competing melting layer models.

Fig. 14.

The 10.7-GHz melting layer optical depth averaged in 1-dBZ bins of basal reflectivity for the FS-a75 snow density and ventilation coefficient combinations of Tables 4 –6. The selected model is shown in bold and the error bars represent the ambiguous region for a wind speed knowledge of 1 m s−1.

Fig. 14.

The 10.7-GHz melting layer optical depth averaged in 1-dBZ bins of basal reflectivity for the FS-a75 snow density and ventilation coefficient combinations of Tables 4 –6. The selected model is shown in bold and the error bars represent the ambiguous region for a wind speed knowledge of 1 m s−1.

Using this technique, the competing dielectric mixing formulas for the melting hydrometeors will be the easiest to resolve. From Table 3, it can be seen that the retrieved wind speed averaged over these profiles varies from 0.25 to 11.89 m s−1 depending on the order in which the dielectric constants of air, ice, and water are mixed using the MG formula. From Tables 4 –6, it is observed that the largest maximum to minimum difference in the retrieved wind speed due to the different parameterizations for the ventilation coefficient for a fixed dielectric formula and snow density is 1.5 m s−1. The largest maximum to minimum difference due to different snow density models for a fixed ventilation coefficient and dielectric formula is 3 m s−1. Thus, in terms of the ability of this technique to constrain the melting layer model, the distinction of the dielectric formula will place the minimum requirements on surface wind speed knowledge, followed by the choice for the snow density model, and then the parameterization for the ventilation coefficient.

9. Summary

A physically based method has been developed to estimate microphysical parameters of the melting layer simultaneously with the rain-rate profile and the underlying surface wind speed using an airborne 10.7-GHz radiometer and a 13.4- and 35.6-GHz radar. The radar backscatter measurements are a function of the drop size distribution in the scattering volume. The frequency dependent scattering properties for the precipitation particles are given by Mie theory. The algorithm employs an iterative inversion technique to estimate the offset parameter, N*0, and the median volume diameter, D0, of the normalized gamma drop size distribution at each radar range gate from the dual-frequency radar backscatter measurements. This requires that the thermodynamic phase of the particles be known. A brightband detection algorithm is used to determine if the profile is convective or stratiform. Without a discernable bright band, the profile is assumed to be convective and the radar profile is divided into a rain and a snow layer, where the transition occurs at the freezing level given by the brightband algorithm. If a bright band exists, the profile is assumed to be stratiform and is divided into a rain layer, melting layer, and snow layer.

To model the backscattering measurements in the melting layer requires a thermodynamic model to determine the meltwater content for a given particle as a function of particle size and height below the freezing level and a dielectric mixing formula to determine the scattering and absorption properties of the particle—assumed to be a mixture of air, ice, and water. The major uncertainties associated with the melting layer model are the parameterization of the ventilation coefficient, the snow density model and the dielectric mixing formula. Several parameterizations for each were evaluated based on their fit to the PR-2 data. It was found that tuning the initial density transition location to 0.75 in the Fabry and Szyrmer (1999) core–shell model provided an excellent fit to the reflectivity peak in the melting layer at both frequencies. This model, when combined with the Magono and Nakamura (1965) snow density model and the Mitra et al. (1990) parameterization of the ventilation coefficient fit both the width and the peak of the PR-2 reflectivity profiles in the melting layer. This melting layer model was used for retrievals in stratiform rain.

The parameterized gamma drop size distribution at each radar range gate was used to determine the extinction and absorption coefficient profile of the precipitation. This was combined with an assumed cloud liquid water and gaseous absorption profiles to determine the 10.7-GHz optical depth and upwelling and downwelling TBs. The radiative transfer algorithm was than inverted to determine the sea surface emissivity, which has a monotonic dependence on wind speed. An empirical forward model from Pandey and Kakar (1982) was used to estimate the near-surface wind speed from the emissivity.

The quality of the surface wind speed estimate is directly related to the quality of the precipitation retrievals. In stratiform rain, the absorption in the melting layer contributes 30% to 50% of the total atmospheric absorption at 10.7 GHz. This means that the quality of the wind retrieval will be driven by the choice of melting layer model. The melting layer model was chosen based on its fit to the PR-2 reflectivity profiles and it was determined that the use of this model produced reasonable wind speed retrievals.

a. Future applications of work presented

An algorithm of this form could be used with the proposed Global Precipitation Measurement Mission. The core satellite of this mission will have a dual frequency radar and multifrequency radiometer. The additional radiometer channels could be used to minimize some of the uncertainties addressed in this study. The largest uncertainty that could be minimized is that in the cloud liquid water content. Retrieval of the integrated cloud liquid water is constrained by the addition of a 37-GHz radiometer channel in light rain and a 19-GHz radiometer channel in moderate rain. Accurate estimates of the precipitation content and underlying surface wind speed in extensive storms would improve the storm-track forecasts by numerical weather prediction models. Sensitivity of errors in the outputs due to errors in the initial inputs for these models double roughly every 2–3 days (Hogan and Rosmond 1991). Current satellite estimates of the wind speed under precipitation from a scatterometer or a radiometer have uncertainties near ±10 m s−1 (Goodberlet et al. 1990; Weissman et al. 2002). The technique demonstrated in this study could be used to retrieve the wind speed under an area of precipitation with an uncertainty of less than 4 m s−1. This assumes that the radar and radiometer are footprint-matched, as was the case in this study. The retrieval accuracy of this technique will most likely degrade for radar–radiometer combinations that have different viewing geometries.

The melting layer model used in this study, which was tuned to fit the PR-2 data at 13.4 and 35.6 GHz, could be used to improve radiometer-only rainfall retrievals. Algorithms that use a Bayesian approach to select precipitation profiles from a modeled TB database will overestimate the amount of precipitation in stratiform rain if a melting layer is not included in the forward model used to populate the database. The errors introduced in the rain retrievals by not properly accounting for the melting layer will be largest in light stratiform rain where the freezing level is close to the surface. Minimizing this source of uncertainty will be especially important for the Global Precipitation Measurement (GPM) mission. Eight constellation satellites, equipped only with multifrequency radiometers, will provide three hour revisit times for locations less than 70° latitude. A majority of the precipitation in midlatitudes is in the form of stratiform rain with low freezing levels. Consistency of the rain rate retrievals between the dual-frequency radar on the core satellite and the radiometer-only retrievals by the constellation satellites will be imperative if the data is to be used in forecast models.

Future analysis of the melting layer model developed in this study should include radiometer observations at multiple frequencies and comparisons in heavier stratiform rain. A field campaign with a multiple frequency airborne radar and radiometer overflying areas of light to moderate stratiform precipitation, preferably with in situ measurements of some of the environmental parameters [i.e., dropsondes, buoys, disdrometers, water vapor radiometers (WVRs)], would be ideal. Even in the current aircraft configuration, the technique developed in this study could be combined with observations of the surface wind speed under stratiform rain as a sensitive and robust means to quantify some of the major uncertainties in the melting layer model.

A long-term ground-based experiment could also be devised to better quantify the major uncertainties associated with brightband modeling. An ideal setup would contain an upward-looking radiometer having 10-, 19-, 21-, and 37-GHz channels and an upward looking dual-frequency radar, possibly Ku–Ka band. A network of disdrometers could be used to estimate the DSD at the surface and a rain gauge grid could be used to constrain the total rainfall estimate. If the forward model for the radar backscatter and the brightness temperature were bound to a physical framework describing the precipitation process, the major uncertainties could be better constrained and the RMS difference between the model and the measurements reduced. The surface DSD and rain rate estimates, as well as surface temperature and humidity measurements, would eliminate some of the unknowns. Long-term statistics could be compiled and used to provide insight into some outstanding questions. For example, can a single model for the dielectric constant of the melting hydrometeors be used in all situations—is it enough to just change the initial snow density, or must the dielectric model be adapted for different situations? What is the role of aggregation and breakup in the melting layer, how does it affect the results, and is it the exception or the rule? Can the bright bands that form in a decaying thunderstorm be modeled in the same manner as those that form in stratiform rain? Long term radar and radiometer statistics of brightband-bearing storms could help answer some of these questions.

Acknowledgments

The authors would like to acknowledge the cooperation of Dr. Eastwood Im and the JPL PR-2 team with the June 2003 field campaign. This work was supported in part by the NASA Earth System Technology Office under Grant NAG5-11609.

REFERENCES

REFERENCES
Atlas
,
D.
,
C.
Srivastava
, and
S.
Sekhon
,
1973
:
Doppler radar characteristics of precipitation at vertical incidence.
Rev. Geophys. Space Phys.
,
11
,
1
35
.
Barthazy
,
E.
,
W.
Heinrich
, and
A.
Waldvogel
,
1998
:
Size distribution of hydrometers through the melting layer.
Atmos. Res.
,
47–48
,
193
208
.
Battaglia
,
A.
,
C.
Kummerow
,
D.
Shin
, and
C.
Williams
,
2003
:
Constraining microwave brightness temperatures by radar brightband observations.
J. Atmos. Oceanic Technol.
,
20
,
856
871
.
Bauer
,
P.
,
J.
Baptista
, and
M.
Iulis
,
1999
:
The effect of the melting layer on the microwave emission of clouds over the ocean.
J. Atmos. Sci.
,
56
,
852
867
.
Bohren
,
C.
, and
L.
Battan
,
1982
:
Radar backscattering of microwaves by spongy ice spheres.
J. Atmos. Sci.
,
39
,
2623
2628
.
Dou
,
X.
,
J.
Testud
,
P.
Amayenc
, and
R.
Black
,
1999
:
The concept of a normalized gamma distribution to descibe raindrop spectra and its use to parameterize rain relations. Preprints, 29th Int. Conf. on Radar Meteorology, Montreal, QC, Canada, Amer. Meteor. Soc., 625–628
.
Fabry
,
F.
, and
W.
Szyrmer
,
1999
:
Modeling of the melting layer. Part II: Electromagnetic.
J. Atmos. Sci.
,
56
,
3593
3600
.
Goodberlet
,
M.
,
C.
Swift
, and
J.
Wilkerson
,
1990
:
Ocean surface wind speed measurements of the Special Sensor Microwave/Imager (SSM/I).
IEEE. Trans. Geosci. Remote Sens.
,
28
,
823
828
.
Grecu
,
M.
, and
E.
Anagnostou
,
2002
:
Use of passive microwave observations in a radar rainfall-profiling algorithm.
J. Appl. Meteor.
,
41
,
702
715
.
Hitschfeld
,
W.
, and
J.
Bordan
,
1954
:
Errors inherent in the radar measurement of rainfall at attenuating wavelengths.
J. Meteor.
,
11
,
58
67
.
Hogan
,
T.
, and
T.
Rosmond
,
1991
:
The description of the Navy Operational Global Atmospheric Prediction System’s spectral forecast model.
Mon. Wea. Rev.
,
119
,
1786
1815
.
Hollinger
,
J.
,
1970
:
Passive microwave measurements of the sea surface.
J. Geophys. Res.
,
75
,
5209
5213
.
Im
,
E.
,
S.
Durden
,
G.
Sadowy
, and
L.
Li
,
2002
:
Rainfall observations by the airborne dual-frequency precipitation radar during CAMEX-4. Proc. 2002 IGARSS, Toronto, ON, Canada, IEEE, 281–283
.
Klaassen
,
W.
,
1988
:
Radar observations and simulation of the melting layer of precipitation.
J. Atmos. Sci.
,
45
,
3741
3753
.
Locatelli
,
J.
, and
P.
Hobbs
,
1974
:
Fall speeds and masses of solid precipitation particles.
J. Geophys. Res.
,
98
,
2757
2765
.
Magono
,
C.
, and
T.
Nakamura
,
1965
:
Aerodynamic studies of falling snowflakes.
J. Meteor. Soc. Japan
,
43
,
139
147
.
Maxwell Garnet
,
J.
,
1904
:
Colours in metal glasses and in metallic films.
Philos. Trans. Roy. Soc. London.
,
A203
,
385
420
.
Meneghini
,
R.
,
H.
Kumagai
,
J.
Wang
,
T.
Iguchi
, and
T.
Kozu
,
1997
:
Microphysical retrievals over stratiform rain using measurements from an airborne dual-wavelength radar-radiometer.
IEEE. Trans. Geosci. Remote Sens.
,
35
,
487
506
.
Meneghini
,
T.
, and
L.
Liao
,
1996
:
Comparisons of cross sections for melting hydrometeors as derived from dielectric mixing formulas and a numerical method.
J. Appl. Meteor.
,
35
,
1658
1670
.
Meneghini
,
T.
, and
L.
Liao
,
2000
:
Effective dielectric constants of mixed-phase hydrometeors.
J. Atmos. Oceanic Technol.
,
17
,
628
640
.
Mitchell
,
D. L.
,
R.
Zhang
, and
R. L.
Pitter
,
1990
:
Mass-dimensional relationships for ice particles and the influence and riming on snowfall rates.
J. Appl. Meteor.
,
29
,
153
163
.
Mitra
,
S.
,
O.
Vohl
,
M.
Ahr
, and
R.
Pruppacher
,
1990
:
A wind tunnel and theoretical study of the melting behavior of atmospheric ice particles. Part IV: Experiment and theory for snow flakes.
J. Atmos. Sci.
,
47
,
584
591
.
Olson
,
W.
,
P.
Bauer
,
N.
Viltard
,
D.
Johnson
,
W.
Tao
,
R.
Meneghini
, and
L.
Liao
,
2001
:
A melting-layer model for passive/active microwave remote sensing applications. Part I: Model formulation and comparison with observations.
J. Appl. Meteor.
,
40
,
1145
1163
.
Pandey
,
P.
, and
R.
Kakar
,
1982
:
An empirical microwave emissivity model for a foam-covered sea.
IEEE J. Oceanic Eng.
,
7
,
135
140
.
Rinehart
,
R.
,
1997
:
Radar for Meteorologists.
Rinehart Publications, 149 pp
.
Rodgers
,
C.
,
2000
:
Inverse Methods for Atmospheric Sounding: Theory and Practice.
World Scientific, 238 pp
.
Rogers
,
R. R.
, and
M. K.
Yau
,
1996
:
A Short Course in Cloud Physics.
3d ed. Butterworth-Heinemann, 290 pp
.
Ruf
,
C.
, and
C.
Principe
,
2003
:
X-band lightweight rainfall radiometer first light. Proc. 2003 IGARSS, Toulouse, France, IEEE, 1701–1703
.
Ruf
,
C.
, and
Coauthors
,
2002
:
Lightweight rainfall radiometer STAR aircraft sensor. Proc. 2002 IGARSS, Toronto, ON, Canada, IEEE, 850–852
.
Rutledge
,
S.
, and
P.
Hobbs
,
1983
:
The mesoscale and microscale structure and organization of clouds and precipitation in midlatitude cyclones. Part VIII: A model for the “seeder-feeder” process in warm-frontal rainbands.
J. Atmos. Sci.
,
40
,
1186
1206
.
Sekhon
,
R.
, and
R.
Srivastava
,
1970
:
Snow size spectra and radar reflectivity.
J. Atmos. Sci.
,
27
,
299
307
.
Stogryn
,
A.
,
1972
:
The emissivity of sea foam at microwave frequencies.
J. Geophys. Res.
,
77
,
1658
1666
.
Szyrmer
,
W.
, and
I.
Zawadzki
,
1999
:
Modeling of the melting layer. Part I: Dynamics and microphysics.
J. Atmos. Sci.
,
56
,
3573
3592
.
Tao
,
W.
, and
J.
Simpson
,
1993
:
Goddard Cumulus Ensemble Model. Part I: Model description.
Terr. Atmos. Oceanic Sci.
,
4
,
35
72
.
Tripoli
,
G. J.
,
1992
:
A nonhydrostatic mesoscale model designed to simulate scale interaction.
Mon. Wea. Rev.
,
120
,
1342
1359
.
Uhlhorn
,
E.
, and
P.
Black
,
2003
:
Verification of remotely sensed sea surface winds in hurricanes.
J. Atmos. Oceanic Technol.
,
20
,
99
116
.
Ulaby
,
F.
,
R. K.
Moore
, and
A. K.
Fung
,
1981
:
Microwave interaction with atmospheric constituents. Microwave Remote Sensing: Active and Passive, Vol. I, Fundamentals and Radiometry, Artech House, 256–343
.
Ulaby
,
F.
,
R. K.
Moore
, and
A. K.
Fung
,
1986
:
Appendix E. Microwave Remote Sensing: Active and Passive, Vol. III, From Theory to Applications, Artech House, 2020–2022
.
Ulbrich
,
C.
,
1983
:
Natural variations in the analytical form of the raindrop size distribution.
J. Climate Appl. Meteor.
,
22
,
1764
1775
.
Viltard
,
N.
,
C.
Kummerow
,
W.
Olson
, and
Y.
Hong
,
2000
:
Combined use of the radar and radiometer of TRMM to estimate the influence of drop size distribution on rain retrievals.
J. Appl. Meteor.
,
39
,
2103
2114
.
Weissman
,
D.
,
M.
Bourassa
, and
J.
Tongue
,
2002
:
Effects of rain rate and wind magnitude on sea winds scatterometer wind speed errors.
J. Atmos. Oceanic Technol.
,
19
,
738
746
.
Willis
,
P.
,
1984
:
Functional fits to some observed drop size distributions and parameterization of rain.
J. Atmos. Sci.
,
41
,
1648
1661
.
Willis
,
P.
, and
A.
Heymsfield
,
1989
:
Structure of the melting layer in mesoscale convective system stratiform precipitation.
J. Atmos. Sci.
,
46
,
2008
2025
.
Wusheng
,
J.
, and
P.
Wang
,
1999
:
Ventilation coefficients for falling ice crystals in the atmosphere at low-intermediate Reynolds numbers.
J. Atmos. Sci.
,
56
,
829
836
.

Footnotes

* Current affiliation: Jet Propulsion Laboratory, Pasadena, California

Corresponding author address: Shannon Brown, Jet Propulsion Laboratory, 4800 Oak Grove Drive, M/S 168-314, Pasadena, CA 91109. Email: shannon.t.brown@jpl.nasa.gov