033112-F EVALUATION OF RADAR TECHNIQUES FOR ASSESSING SNOWCOVER CONDITIONS AND THEIR EFFECT ON THE DETECTION OF HARD TARGETS FINAL REPORT K. Sarabandi and F.T. Ulaby Office of Naval Research Code 321 800 North Quincy Street Arlington, VA 22217 Contract: N00014-95-1-0736 June 1996 33112-1-F = RL-2454

033112-F EVALUATION OF RADAR TECHNIQUES FOR ASSESSING SNOWCOVER CONDITIONS AND THEIR EFFECT ON THE DETECTION OF HARD TARGETS Final Report S.O. Code 321 Office of Naval Research Contract N00014-95-1-0736 June 1996 By: K. Sarabandi and F.T. Ulaby THE VIEW, OPINIONS, AND/OR FINDINGS CONTAINED IN THIS REPORT ARE THOSE OF THE AUTHORS AND SHOULD NOT BE CONSTRUED AS AN OFFICIAL DEPARTMENT OF THE NAVY POSITION, POLICY, OR DECISION, UNLESS SO DESIGNATED BY OTHER DOCUMENTATION.

1 INTRODUCTION The overall objective of this investigation was to determine how well a down-looking radar system can detect stationary hard targets against a snow-covered terrain background. The presence of snow cover and its condition is of importance for Navy as it may pose trafficability problems for vehicles and supplies, and therefore knowledge of the spatial distribution of snow depth and related snow parameters over an area of interest could prove very useful in deployment exercises. Over the past two years, extensive theoretical, and experimental efforts have been devoted to the understanding of the polarimetric response of snow-covered terrain at microwave frequencies and to assess the applicability and accuracy of the existing models for retrieval of snowpack parameters. Our experimental efforts were directed at: (1) conducting polarimetric backscatter measurements for snowpacks over a wide range of incidence angles and at three frequencies 1.25 GHz, 5.3 GHz, and 9.5 GHz, (2) development of sensors for "ground-truth" measurements and collecting accurate snow parameter data. The experimental data were used to assess the accuracy of the existing theoretical and empirical models. We have also been able to develop a new approach to model dense random media such as snow. This approach is based on combining the experimental and theoretical scattering behavior of the medium of interest where the fundamental quantities of the Radiative Transfer model (phase and extinction matrices) are retrieved experimentally instead of analytically. The outcome of our research activities are presented in form of two journal papers which are given in Appendices A and B. These papers were submitted to the IEEE journal Transactions on Geoscience and Remote Sensing. The first paper entitled "Radar measurements of snow: Experiment and analysis" provides an extensive literature survey of existing theoretical and empirical models for snowpacks and includes the details of the radar experiments conducted for snow-covered ground surfaces. In this study it was found that both conventional and dense-medium radiative transfer models fail to adequately explain the experimental results. It was also shown that polarimetric scattering and extinction quantities intrinsic to snow media can be retrieved accurately which can then be used to predict the scattering behavior of the sowpack. A semi 1

empirical model (reference [8] Appendix A) was also examined against the experimental data and was found that the semi-empirical model was able to predict the snow wetness reasonably well. The second paper entitled "A hybrid experimental/theoretical model for dense random media" outlines a new approach to model the scattering and propagation behavior of electromagnetic waves in complex dense random media such as snow. This paper includes polarimetric backscatter measurements made at Ku-band on layers of a dense medium under very carefully controlled circumstances which are used to: (i) evaluate the degree to which the experimental observations are predicted by theoretical, particle-based, random media models; and (ii) test the proposed hybrid model. The hybrid model employs the mathematical machinery of the first-order vector radiative transfer. In this study very simple forms of the phase matrix components involved in the first-order formulation have been derived using Rayleigh theory. The major conclusions of this study are that: (i) the hybrid model is an adequate description of the dense medium scattering behavior, (ii) conventional radiative transfer (RT) appears to give a better estimate of the observed radar response than the dense medium RT; and finally, (iii) the phase function of the effective volume scattering element of the medium, obtained via the hybrid model, suggests a larger effective scatterer than the physical ones. 2

Appendix A Radar Measurements of Snow: Experiment and Analysis 3

Radar Measurements of Snow: Experiment and Analysis John R. Kendra, Kamal Sarabandi, and Fawwaz T. Ulaby Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Abstract This paper considers two specific types of experiments conducted to improve our understanding of radar backscatter from snow-covered ground surfaces. The first experiment involves radar backscatter measurements at C- and X-band of artificial snow of varying depths. The relatively simple target characteristics, combined with an exhaustive ground truth effort make the results of this experiment especially amenable to comparison with predictions based on theoretical methods for modeling volume scattering media. It is shown that both conventional and dense-medium radiative transfer models fail to adequately explain the observed results. A direct polarimetric inversion approach is described by which the characteristics of the snow medium are extracted from the measured data. The second type of experiment examined in this study involves diurnal backscatter measurements that were made contemporaneously with detailed measurements of the snow-wetness depth profiles of the observed scene. These data are used to evaluate the capability of a recently proposed algorithm for snow wetness retrieval from polarimetric SAR measurements, which has hithertofore been applied only to data from very complex and extended mountainous terrains. 1

1 Introduction A great deal of experimental and theoretical work has been done pertaining to the radar response of snow or,by extension, dense random media in the theoretical realm, for which snow is perhaps the most representative natural example. The aim of all of this work, and, what is true of any remote sensing research, is to develop the capability to characterize in some way, to a greater or lesser degree, the remote target from the sensor response(s) alone. We may identify three levels of such characterization, which we list in order of increasing power of characterization but decreasing level of reliability. * Level 1: Empirical models are used to infer or predict information about the target characteristics. Recent examples include radar and radiometer algorithms for discriminating wet snow covered terrain from other types of terrain [1], a hybrid empirical/theoretical approach for estimating radar clutter at millimeter wave due to certain types of terrain [2], and an approach for classifying snow cover states (dry/wet/refrozen) [3]. * Level 2: Classification is performed based on the apparent scattering mechanisms inherent in the target as identified through analysis of the character of the target Mueller matrix [4, 5]. Such techniques have reportedly been successful in demarcating areas consisting of primarily urban targets, slightly rough targets like oceans and lava flows, and parks and vegetated areas. It has also been shown that such a method may allow discrimination between relatively younger and older lava flows [6]. These techniques have also lately been suggested as a means for determining wetness levels in snow [7, 8]. * Level 3: At this level, tools from the previous two may be used, but the central characteristic is the use of a theoretical model which is assumed to generate reasonably high fidelity predictions of sensor responses for a given set of physical parameters which is assumed to constitute an accurate description of the target. An example of this type of approach is a neural-net-driven inversion algorithm intended to allow retrieval of snow parameters, from radar and radiometer sensor responses, and which is trained using a dense medium radiative transfer model [9]. It is obvious from the descriptions that, where greater understanding of the electromagnetic interaction with a material is present, the potential for information retrieval through remotely sensed data is greater. The critical issue becomes then testing the validity of theoretical models through careful experiments. 2

In the present study, an attempt is made to address this issue with respect to the radar response of snow. There exists already numerous experimental studies of snow in the literature, both at microwave frequencies [10, 11, 12, 13, 14], and at millimeter-wave frequencies [15, 16]. Using such experimental data for the purpose of evaluation of models is very difficult because of the need to carefully characterize the target. Even with precise characterization, substantial obstacles remain. A complex target can be described down to the finest details, but this still leaves the problem of correctly modeling the behavior of all these features, electromagnetically speaking, and all of their interactions with each other. Consider a typical target of snow-covered ground. Parameters which must be potentially considered are: * Parameters associated with the top surface of the snow (rough surface) * The snow itself: density, particle size distribution; vertical distributions of these properties within the snowpack. * Snow wetness, when present, may well be a very complex function of time and depth. Examples of this are presented in later sections of this paper. * The ground beneath: its dielectric constant, roughness parameters, local slope. For experiments at millimeter-wave frequencies, the number of parameters may be somewhat reduced since it is generally accepted, and has been demonstrated by careful experiments [17], that the total penetration depth is relatively shallow. Thus, at mm-wave frequencies, the snow may be treated as a half-space. The price for this, however, is that sensitivity to snow microstructure and other smaller scale features is heightened, for which accurate characterization may be very difficult. Microwave frequencies offer the most potential for the retrieval of gross snow properties such as depth or water equivalent, parameters which are especially important for hydrological applications. For microwave frequencies, generally most of the parameters mentioned above need to be considered. Thus it is in practice difficult to address theoretical predictions in an unambiguous fashion with field experiments. The present study describes radar backscatter experiments on snow at C- and X-band. In these experiments, the use of artificially produced snow allowed an unusually high degree of confidence in the exact character of the target snowpack. Two different kinds of experiments were conducted: experiments of the angular response of dry snow at various depths, and experiments at a single angle over a partial diurnal cycle during which complete profiles of water content as a function of depth and time were recorded, using the 3

Snow Probe [18]. A complete description of the experiment is presented in Section 2. In Section 3 the results of the backscatter experiments are presented and analyzed with respect to certain models. In Section 4, the results of this diurnal experiment and another previous diurnal experiment are presented and discussed with respect to a proposed inversion algorithm for snow wetness. 2 Experiment Description This section describes scatterometer experiments which were performed at C-, and X-band (5.3, and 9.5 GHz) on snowpacks comprised of artificial snow. The experimental site was the Mt. Brighton Ski Area, in Brighton, Michigan, during February and March of 1993. The radars used were truck-mounted, fully polarimetric network-analyzer-based systems, using horn antennas. Detailed descriptions of the systems are available elsewhere [19, 20]. The target was an area of ground covering approximately 18m x 30m. On the average, sixty independent spatial measurements were taken with each angular measurement to reduce the variance due to fading in the estimate of the mean backscatter. Additional frequency averaging was available from the bandwidths used in the respective channels, 400, and 500 MHz in the C-, and X-band channels respectively. Calibration was performed using a fourteen-inch sphere and a differential Mueller matrix algorithm suitable for measurements of distributed targets [21]. Two types of experiments were performed. The first consisted of backscatter measurements at incidence angles of 20~ through 60~ on bare ground, which were then repeated for three progressively deeper layers of artificial snow over the ground. This sequence was performed over a six day period, from February 24th to March 1. A major goal of the experiment was to simplify, insofar as was possible, the nature of the target to allow as unambiguous comparison with theory as possible. One goal we had was insuring that there would be no wetness present in the target snowpacks, so that we would not have to deal with this complicating factor in our interpretation. We therefore conducted all of the experiments during the night. As may be seen from the daily temperatures plotted in Figure 1, the temperature during the data-takes averaged -15~ C and was always below -10~ C. Snow probe measurements verified that there was, to the accuracy of the measurement, essentially no liquid water present in the snowpack. An additional fortuitous circumstance, which may also be observed in Figure 1, was that the temperture was very low in the week preceeding the experiment, leading 4

to a generally frozen environment which was unrelieved during the entire experiment as even the temperature highs did not go above the freezing mark until the day following the final experiment. The bare ground target was prepared by using earth moving equipment to remove all of the existing snow, natural and artificial, and then to scrape away any grass or other vegetation on the target plot. The dielectric constant of the ground was measured using a microwave field-portable dielectric probe [22] operating at C-band. The average dielectric was found, from twenty-five separate measurements to be 4.6 ~ 1.1 (the imaginary part is not estimated accurately through this technique). The bare ground roughness was measured using a laser profilometer. Seven 95-cm linear transects were measured; the step size was 3mm. The average rms height of the surface was found to be, s = 0.32 ~ 0.08 cm; and for the correlation length, I = 2.09 ~ 1.6 cm. The artifical snow was added in three progressively deeper layers. No chemical additives were used in the production of the snow, as is commonly done to increase production, which might have caused the dielectric constant of the particles to depart from that of water ice. The snow, once produced, was spread into a uniform layer using a grooming machine (piston bully) with extra wide treads to reduce the pressure exerted by the weight of the machine on the snowpack. A consequence of this process, was that a characteristic groove pattern was imparted to the top of the snow, having a period of 3.2 cm and a peak-to-peak amplitude of 1.3 cm, corresponding to an rms height of 0.449 cm. The correlation length of such a surface is dependent on the direction with which it is computed relative to the grooves. The three different depths of snow deposited were 20,60, and 102 cm, with the standard deviations shown in Table 1. The nature of artificial snow is that it consists of very small particles and has a very high density. The average density of the snow was 0.48 g/cm3, and the average particle diameter was found to be 0.27 mm. This combination of particle size and density are at the extreme ends of what would be found in nature, but not to the point of being unreasonable. We note, for example that the density is very comparable to that reported for Alpine snowpacks in [7]; neither are these particle sizes outside of reported ranges [23, 17]. All of the target physical properties for the dry snow backscatter experiments are summarized in Table 1. It deserves to be emphasized that there were a number of elements present in this experimental effort which make it unusually well-and simply-characterized. The snow machines produced very uniform "snow" over the course of the experiment, and this was recorded and photographed. From Figure 2, it is evident that the snowpack was in fact composed of discrete particles having a high degree of sphericity. Furthermore, 5

Table 1: Dry Snow Physical Properties Ground Snow Volume Snow Surface Snow Depths s=0.32cm d= 0.27 ~0.11 mm s= 0.45cm di =20~6cm 1=2.09 E = 1.9 - jO.015 e' = 1.97 d2= 60 10 cm - = 4.7 (C-band) Ps = 0.48 g/cm3 Period: 0.32 cm 43 = 102 ~9 cm ____ __________________ Amplitude (p-p): 1.27cm s = rms height, I = surface correlation length, d = mean particle size there was little or no variation from this snow composition as a function of vertical position within the snowpack. This last point is owing not only to the uniform way in which the particles were produced but also to the unbroken cold conditions which persisted throughout this experimental phase and therefore prevented any metamorphic activity associated with melting and re-freezing within the earliest deposited snow. Finally, the action of the grooming machine on the top surface of the snow guaranteed that virtually identical roughness was present for all target realizations. That the surface was not isotropic is not a desirable characteristic, but, as we will show, the surface, having such a small dielectric constant, does not make an especially large contribution to the backscatter anyway. What is most important is that it cannot be invoked to explain differences in the responses of the various target configurations. The second type of experiment performed was of a partial diurnal cycle. This experiment took place on March 6, five days after the last experiment in the dry snowpack sequence. In this intervening period, the daily high temperatures had been consistantly above freezing, getting as high as 9~ C. As a result the snow pack had experienced considerable metamorphosis. In addition, there had been deposited about 15 cm of new snow. Thus, the snowpack had begun to resemble a natural one, complex and difficult to characterize, and mainly unsuitable for comparison with rigorous modeling approaches. However, since copious data on wetness, as a function of both vertical position and time were recorded, along with the polarimetric backscatter at 40~ incidence, it provides an excellent opportunity for direct evaluation of a recent algorithm [7, 8] for inversion of snow wetness levels from target Mueller matrices. A second data set which was recorded on March 1, 1993, in Cadillac, Michigan, also at 40~ incidence and with complete wetness data, will also be considered. 6

3 Dry Snow Backscatter Results & Discussion 3.1 Bare Ground To model the backscatter from the dry (artificial) snow, it is first necessary to account for the contributions of the soil surface below the snowpack as well as the top surface of the snowpack. In [20] it is shown, as the result of an extensive experimental effort, that classical models such as the small perturbation model and the physical optics and geometric optics solutions of the Kirchhoff approximation do not accurately predict the co- and cross-polarized backscatter from bare soil surfaces. The result of that study is a simple semiempirical model (and inversion algorithm) for co-polarized and cross-polarized backscattering coefficients, which uses two parameters: ks, the surface roughness, and m, the soil moisture content. The value of m, for the soil is found to be 0.085 based on our measured soil dielectric at C-band in conjunction with empirical expressions relating these two quantities which are provided in [24]. Using this value, along with the soil rms height s (Table 1), we may compute predictions of the soil backscatter using the semiempirical model from [20]. Using these exact values, the X-band response is somewhat underestimated; however, a small adjustment of the rms height from the measured value of 0.32 cm to 0.37 cm, gives close simultaneous agreement for both frequencies for both co- and cross-polarized responses, as shown in Figure 3. These values, m, = 0.085 and s = 0.37 cm will be used to compute the backscatter for the snow-covered ground. The semi-empirical surface scattering model from [20] will also be used to compute the backscatter contribution from the top surface of the snow. 3.2 Dry Snow Backscatter The results from the angular measurements made on three different depths of artificial snow are shown in Figure 4. Also shown, to better illustrate the snow volume contribution, are simulations of the surface backscatter from (i) the snow-covered ground alone and (ii) the snow-covered ground plus the contribution from the top snow surface. These are shown separately to illustrate the relative importance of the top-surface term. One obvious characteristic of the data is the relatively small dynamic range of the co-polarized responses corresponding to different depths. Indeed the co-polarized response for the shallowest (20 cm) depth can be attributed to the effects of the two rough surfaces alone. Still, there is an upward trend apparent, as the 60 cm and 102 cm depths are on the order of 3 dB higher than the 20 cm case. There appears to be saturation 7

occurring, as the responses from the two deeper layers appear to be essentially the same, within the apparent noise level of the measurements. The level of noise within the measurements appears to be appreciable; this is attributed to the relatively difficult conditions (cold, darkness, finite boundaries of the target area) under which they were made. It is worth comparing these results qualitatively with what has been experimentally observed before for snow in this frequency range. One issue which has been debated in the literature is that of whether snow is observable at all at X-band. Several investigators have claimed that it is not, based on their measurement campaigns [12, 13], whereas, in another case, [11], evidence to the contrary is reported. A critical element which must be taken into account in considering this question is the magnitude of the ground scattering. For example, in [12] and [13], the reported rms roughness of the bare ground was 2 cm. Accordingly, the average bare ground scattering level is as high or higher than the maximum levels measured in our present experiments. In the dry snow-pile experiments of [11], the bare ground characteristics are not reported, but o~ levels for the shallowest snow depths (from Fig. 4 in [11]), are suggestive of a relatively smooth surface, probably similar to that of the present study. The dynamic range of the cross-pol is seen to be considerably greater than the co-polarized, and is unquestionably increased with snow depth, exhibiting an increase on the order of 7 dB. Here too is observed an apparent saturation as in the case of the co-polarized, for the two deepest layers. 3.2.1 Comparison with Discrete-Particle-Based Theories The well-defined nature of the artificial snowpack allows direct comparison with theories which are based on discrete particles. Since intuition suggests that for particles of this size (d = 0.27 mm) the scattering will be small, we will start by considering an independent scatterer formulation, since it is simpler. Since it has been found both experimentally [25] and theoretically [26] that correlated dense medium scattering is less than independent scattering at low frequencies, the solution will represent an upper bound on what dense medium radiative transfer (DMRT) solution techniques may predict. For a comparison with theory, the complex dielectric constant of the ice particles ice, must also be known. For the real part, e',, we use 3.15. The imaginary part is computed using the following formula [27] which 8

compares very favorably with published data and also accounts for temperature dependance: eC = 57.34(- +2.48 x 10 /f)exp(3.62 x 10-2T) (1) where f is in Hz and T is in degrees Kelvin. For particles of this size, the scattering albedo, 0o, as computed using the Mie solution, is only 5.2 x 10-3 for C-band and 2.3 x 10-2 for X-band. The scattering albedo is the ratio of the scattering cross-section to the total extinction cross-section (as/ a). That it is so small in this case indicates that volume scattering may be treated as a perturbation on the reduced (by extinction) coherent wave in the medium. Thus a first-order radiative transfer solution is appropriate for the solution of a layer of these particles. If this solution is computed for a layer of particles (having a smooth surface) over a smooth dielectric half-space having the same dielectric constant as the soil in the present study, an estimate of the contribution of the snow volume is obtained. When we performed this calculation, as a function of depth and angle for the two frequencies, we found that the maximum co-polarized a~ produced was -48 dB at C-band and -32 dB at X-band. The contribution to the total scattering represented by these levels does not appear as a visible increase relative to the curve in Figure 4 which corresponds to the backscattering contributions of the two rough surfaces, top and bottom, alone. The cross-polarized response corresponding the first-order solution for spherical particles is identically zero. Use of a more sophisticated radiative transfer solution [28] which uses the discrete ordinate method and so accounts for all orders of scattering gives, not surprisingly, essentially the identical results for the co-polarized cases; the cross-pol estimate is on the order of -80 dB for C-band and -70 dB for X-band. It is obvious that the behavior of the target cannot be explained in terms of the particles of which the snowpack was observed to be comprised. There has appeared in the literature recently work which considers "sticky" particles [29], that is, particles which come together to form an aggregate particle, effectively much larger than the individual particles. Obviously, this could have a profound effect upon the scattering behavior of a medium, recalling that in the Rayleigh regime, scattering increases (on a per unit volume basis) as r3. If the particle size is treated as a free parameter in the conventional RT formulation, a value can be obtained which gives the optimal agreement with the measured results. The results of such a process are shown in Figure 5. As shown, reasonable agreement-for the co-polarized responses-can be obtained if an effective particle diameter of 1.7 mm (more than six times larger than the measured average diameter of 0.27 mm) is 9

used for C-band and a diameter of 0.9 mm (more than three times larger than the measured value) is used at X-band. Even in this case however, the RT formulation fails completely to predict the substantial crosspolarized response observed in the measurements. Faced with the failure of existing particle-based theories to explain the experimental observations, we turn to the question of what useful information can be be extracted from this data. A record of the snowpack target that was measured adds, in itself, very little to the study and practice of remote sensing of snow. There is a low probability that the same combination of physical features, including the characteristics on top of the snowpack and beneath it, will be duplicated elsewhere. A potentially very useful result may be obtained, however, if the intrinsic quantities which specify the extinction and scattering characteristics of this snow material can be retrieved. This concept is illustrated in the following section. 3.2.2 Direct Characterization of the Snow Medium This section describes an approach for retrieving extinction and scattering parameters for the snow used in this study. A major assumption which is made is that dense media scattering behavior can be described by a first-order radiative transfer formulation. The familiar four terms which results from such a formulation are depicted in Figure 6. We have intentionally represented the scattering elements as clusters to underscore the point that we are considering "effective" particles in this treatment, which may comprise correlated groups of individual physical particles and/or multiple scattering effects. The validity of such an assumption, that is, that dense media scattering can be understood in terms of a first-order radiative transfer model, has been demonstrated in [30]. A second assumption which will be made, which will greatly simplify the analysis, is that only the direct backscatter term (term (A) in Fig. 6) in the volume scattering formulation is important. This assumption is justified due to the relatively small reflectivity associated with a reflection at the snow / soil interface, which is less than -12 dB for either polarization at any of the angles examined. In this case the first-order solution of the radiative transfer equation reduces to the form, Lm = Lis + i 21 Pbs + Lg T12 (2) 10

where, Lis is the Mueller matrix corresponding to the top surface scattering, Tpq is the surface intensity transmissivity matrix [31][p. 145] for transmission from p to q, where 1 corresponds to the air and 2 to the snow medium, g' is the cosine of the refracted angle, Pbs is the unknown backscatter component of the phase matrix for the effective volume scattering element, Lg.d corresponds to the rough surface below, and, 1 - exp(-2Ked/g') 2K1e/g L2 = exp(2Kced/g'), where Y, is the extinction and d is the layer depth. This model has in the past been applied in a scalar sense to both snow [11] and vegetation [32]. In order to apply (2) in a vector sense, it is necessary to construct Mueller matrices corresponding to the top and bottom rough surfaces. In the preceeding analysis, we had used a semi-empirical model for rough surface scattering because of its greater reliability. The model as it has been presented is a scalar one; in [33], however, Oh extends this approach to include an empirical model for a, the degree of correlation between S, and Shh. This parameter governs [34] the width of the probability density function which describes the distribution 4hh - (v,, the phase difference of the co-polarized complex scattering amplitudes. Oh's formula is given as: a = 1 [ - 0.2(sin0)A(ks,)] (cos 0)B(kso) (3) where, A(ks,ro) = (16.5ro +5.6)exp[-41.6ksF2] (4) B(ks,1ro) = 8.1roksexp[-1.8ks] (5) where ro is the Fresnel reflectivity at nadir, k is the free-space propagation constant, and s is the surface rms height. The importance of knowing a for the present purpose is that it is expressible in terms of the elements of the Mueller matrix, and thus can aid us in constructing these matrices for the ground and the snow surface. 11

The quantity a is formally given by [34] as, [- = 1 MJ 12 (6) [ (M33 + M44)2 + (M34-M43 )2]/2 (6) where Mm,, are elements of the modified Mueller matrix, L4: < ISvl2 > < IShh12 > < Re(S*hSVV) > < ISvl2 > < IShh|2 > < Re(S^hShv) > 4nLm < 2Re(SvvSv) > < 2Re(ShSh)) > < Re(SvvSwh + SvhSv) > < (SS)> 2Im(S ) >(SvhSSh) > < -Im(SvS^ -SvhS*Y) > < -Im(S*Svv) > < -Im(SShh) > (7) < -Im(SvvS;h -SvhS*v) > < Re(SvvSh-S vhS;v) > From Eq. (6), if we make the assumption that M43, which corresponds (for backscatter) to the term - < Im(SwS^h >, is much smaller than < Re(SvvSh) >, then < SvvSh^ > ~avM11M22 (8) From comparison with the bare ground backscatter data, the selection of the negative root is indicated, (following the FSA convention), which is also in agreement with polarimetric predictions for rough surfaces [4] in which it has been found that the statistical phase difference between Svv and Shh is near 180~ (0~ in the BSA convention). For the bare ground under the snow layer, a was found, using (3) to be essentially unity for both C-band and X-band for all angles between 20~ and 60~ degrees. For the top snow surface, a varied from unity at 200 to 0.89 at 600 for both frequencies. A further assumption that was made in constructing the Mueller matrices for the top and bottom surfaces-and one which is reasonable given the above discussion regarding the co-polarized phase difference for rough surfaces-is that the elements M34 and M43 corresponding to ~Im < SVVS^h > are zero. 12

The polarimetric character of the two rough interfaces, above and below the snow layer are thus characterized by a target Mueller matrix having the following form: gvv gvh 0 0 gvh ghh 0 0 Lsurf = (9) 0 0 -aO/gvvghh - gh 0.0 0 0 -/gvvghh + gvh where, 0~ gpq - 4cos0, (10) with oGq as given by the semi-empirical model of [20]) and On the angle of incidence onto the surface in the appropriate medium (air for the top surface, snow for the underlying soil surface). The correlation coefficient a is as calculated in Eq. (3). The parameters to be specified through an optimization of the measured polarimatric data with the model given in (2) are the scalar extinction and the backscatter component of the phase matrix 2b,. Since the snow is assumed to be an isotropic medium we use the form: P1 P2 0 0 P2 P1 0 0 Pbs = (11) 0 0 P3-P2 -P4 0 0 P4 P3 + P2 which requires the volume backscattering for the two co-polarized channels to be identical. The nature of the other elements is set by considerations of reciprocity (ISvh12 = ISh v2) and zero correlation between the coand cross-polarized complex scattering amplitudes. With the unknown scalar extinction there are therefore five unknown quantities to be determined in order to characterize the snow polarimetrically. The estimation of these parameters is done using an optimization package suitable for non-linear problems [35]. The optimization is performed using only the 40~ data, considering all three depths, for both C-band and X-band, and then the results applied to generate predictions for the other cases (depths and angles) not used in the 13

Table 2: Artificial snow: Empirical parameters for RT-model Parameter C-band X-band e (Np/m) 0.513 1.28 Pi (m-1) 0.316 x 10-2 0.128 x 10 P2 (m-') 0.330 x 10-3 0.740 x 10-3 P3 (m-1) -0.169 x 10-3 -0.976 x 10-3 P4 (m- ) 0.531 x 10-3 0.220 x 10-2 Table 3: Estimates of the extinction of the artificial snow (a, Np/m) Source C-band X-band "Meas." 0.513 1.28 EFA 0.020 0.047 QCA-CP 0.042 0.096 optimization process. In Figures 7 thru 9, we show the measured data along with predictions generated using the estimated parameters Ke and P1,...,P4 in (2). The angular variation predicted by (2) (in which the angular variation is contained in the transmissivity matrices) appears reasonable, as does the behavior with respect to depth at angles which were not involved in the optimization process. Thparameters Ke and Pi,...,P4 for both C- and X-band are summarized in Table 2. A comparison of the values for extinction obtained by this empirical and two theoretical methods is given in Table 3. One of the two methods, EFA or the Effective Field Approximation [36], is associated with conventional radiative transfer (CRT) and the other, QCA-CP or, Quasi-Crystalline Approximation-Coherent Potential [36], is used in dense medium radiative transfer (DMRT). It is seen from Table 3, as for the scattering computations mentioned earlier, the predicted effects of the snowpack are practically negligible. Since the analysis was a complete polarimetric one, we can examine the results with respect to the two remaining (given the two co-polarized responses and the cross-polarized response already examined) independent quantities associated with a measured Mueller matrix in the backscatter direction. These quantities are the Degree of Correlation, cc, which was introduced earlier, and another quantity which pertains to the position of the maximum of the probability density function [34] describing the distribution of the phase difference between Svv and Sh^. This quantity is known as the Co-polarized Phase Difference and is denoted by the symbol 4. Figure 10 compares the measured values of these quantities at C-band and the corresponding estimates 14

generated throught the use of the RT model with parameters given in Table 2, for the case of the 60 cm artificial snow layer. The agreement is not exceptionally good. This is attributed to (i) the need to specify the polarimetric character of both rough surface (above and below the snow) and (ii) the generally noisy character of the data. 3.3 Discussion of Dry Snow Results The most important result from the dry snow experiments was the observation that the most widely used models describing the electromagnetic scattering characteristics of dense media failed to give reasonable predictions in comparison to measured data. While one experiment is not sufficient for evaluating the valdity of a theory, the present experiment carries significant weight because of the relatively simple and well-characterized nature of the target snowpack. In the subsequent analysis, a simple method for retrieving electromagnetic parameters intrinsic to the snow medium itself was described. It is noteworthy-in the face of the very large body of research which has been done on snow and pertaining directly to snow-that this marks the first time an attempt has been made to characterize the effect of a snow medium considered in isolation from other effects, in a polarimetric way. The only other comparable example is the snow-pile experiments and subsequent scalar analysis described in [ 11], in which the radar response of a snow-pile of varying depths was collected. That effort was hampered by a number of departures from ideal experimental circumstances. For example, only one spatial sample was available, the snow was artificially heaped up by mechanical means, and efforts to characterize the upper or lower surfaces or even the general uniformity of the snow pile were absent. Of course, for direct characterization of snow, it is not feasible to use the approach employed in this experiment. A practical scheme should borrow from the spirit of the snow-pile experiments alluded to above, in terms of having a target of a fairly manageable size-except that a turntable could be employed to allow realization of independent samples and steps taken to control the character of the surfaces above and below. The material representing the underlying half-space could be chosen to best facillitate the retrieval of information. Feasibility studies on this technique, using not snow but stable materials like sand and gravel, have in fact been carried out by the author and the results, including an analysis of the validity of a vector radiative transfer model with empirically-derived parameters for describing very dense media may be found in [30]. 15

4 Diurnal Results and Discussion This section present the results from measurements of partial diurnal cycles which were collected on two separate occasions, as described in the Section 2. The snow probe was used on both occasions to record vertical profiles of the liquid water content as a function of time. The two snowpacks examined had very different physical descriptions as will be shown. This allows some insight into the generation and spatial behavior of liquid water in a snowpack. We will use the diurnal results and associated wetness data to evaluate a particular algorithm which has been recently proposed for the retrieval of snow liquid water content. 4.1 Brighton Diurnal Results In the Brighton diurnal experiment, the measurements were made continuously from 10 a.m. until 7 p.m.. Overnite temperatures prior to the experiment were well below freezing; the temperatures during the day ranged from -61C at 8 a.m. to more than 6VC at 3 p.m. It was a very sunny day, and heavy melting was evident. By 6 p.m. the temperature dropped abruptly below freezing again. The target snowpack was 0.88 -m deep, with about 15-cm of relatively fresh snow on top, however, the bottom 16 cm consisted of solid ice. The density of the snowpack was about 0.25 g/cm3 at the top, increasing linearly to 0.45 by 30 cm into the pack, where it remained essentially constant to the bottom of the snowpack. The co- and cross-polarized results for C- and X-band are shown in Figure 11. Both frequencies show a significant reduction in the backscatter at midday, very typical of a diurnal response, in this case 8-10 dB for C-band, and - 14 dB for X-band. By the end of the measurements however, the radar response at both frequencies appears to be headed back up towards the original morning levels. The results of the snow probe measurements made concurrently with the radar measurements are shown in Figure 12. The wetness measurements begin at 18-cm, due to the presence of the ice layer below this, and are made at roughly 5-cm intervals. The temporal spacing between measurements of the vertical profile was ~ 1 hour. What is most striking about the wetness map is the presence of very significant wetness levels in the lower 35-40 cm of the snowpack, even at the earliest point measured in the morning, while the top surface is completely dry. The wetness level of the top surface remains fairly moderate through the day, staying below 5% except, curiously, at the very end of the day, just before the temperature fell very swiftly below freezing again. It would appear that the top of the snowpack was freely draining throughout the day as the wetness levels towards the bottom 16

were observed to increase to maximal snow wetness levels, >12%. 4.2 Cadillac Diurnal Results In the Cadillac diurnal experiment, the snowpack was only 22 cm deep. The measurements were made from 10 a.m. to about 4 p.m. As in the Brighton case, the ovemite temperatures had been sub-freezing but the experiment day itself became very warm, reaching 120C for a high, and very sunny. Extensive melting was clearly evident around the entire experimental area. The co- and cross-polarized results are shown in Figure 13, and the associated wetness map in Figure 14. The apparent higher frequency of collection of radar measurements relative to the Brighton data pertains to the number of measurements which were combined as a single data point in each case. These results are seen to have a very different character than that observed in the Brighton case. An important feature of the Cadillac snowpack was the presence, at least during the initial few hours of the experiment, of a very prominant ice lens, starting approximately 2.5 cm below the top surface and having a thickness of about 2 cm. A result which is attributable to both this feature and the relatively much warmer temperatures which occurred compared to the Brighton case, is the presence of very high wetness levels in the uppermost levels of the snowpack. Though the high wetness levels are not strictly confined to the 2.5 cm above the ice lens, it is apparent that the ice lens, particularly early on and to a lesser extent as it became softer and more permeable, impeded the drainage of the liquid water through the snowpack. Indeed there appears to be some evidence from an examination of Figure 14 that the ice lens rapidly evolved, as it softened, into a region supporting very high wetness levels. One implication of this situation (high wetness levels near the surface) evident in the backscatter (Fig. 13) is a "hump", or at least a temporary departure from the downward trend of o~, for both co- and crosspolarized channels. This feature, occurring between the 11:30 a.m. and about 2 p.m., is apparently tied to the increase in surface scattering, due to a higher dielectric contrast, balancing the reduction in volume scattering. 4.3 Evaluation of Wetness Retrieval Algorithm The availability of detailed snow wetness information along with polarimetric backscatter data allows comparison with an algorithm which has recently been developed [7, 8], for the retrieval of snow liquid water estimates from polarimetric data at C-band. The ability to detect, quantitatively, from remotely-sensed data, 17

a parameter such as liquid-water content in snow would constitute a very important achievement and a major step forward in supplying hydrologists and other earth scientists with data products critical for various applications. It is all the more impressive given the difficulties associated even with direct measurements of liquid water content. Given the potential value of such an algorithm, it is highly desirable that its validity be supported by careful ground-based measurement efforts. The present diurnal experiment appears to be an excellent candidate for this evaluation. The algorithm was originally motivated by an AIRSAR data set which was collected over glaciers in the Oztal Alps, Austria. Prior to these measurements, physical characteristics such as snow depth, density, wetness and surface roughness were measured. The final version of the algorithm, as appears in [8] was developed in conjunction with data from a recent (Apr. 1994) SIR-C/X-SAR mission, taken over Mammoth Mountain on the eastern slope of the Sierra Nevadas. A similar ground truth campaign was employed in this effort as well. The algorithm conceives of the radar response as an incoherent addition of volume and surface responses. No attempt is made to explicitly model the magnitude of the volume contribution to the scattering. Instead, it is shown how under certain assumptions, the relationships between certain quantities in the polarimetric volume radar response may be explicitly constrained. Specifically, the volume, since wet snow is assumed, is considered a half-space. A first-order volume scattering mechanism is assumed which leads to the result that the ratio of the co-polarized volume responses is equal to the square of the ratio of the Fresnel transmissivities corresponding to each polarization respectively. Similar ratios are constructed between the co-polarized volume responses and a term associated with the correlation between the two co-polarized channels: Ovvhh = (Re[SvShh]). (12) The correlation is assumed to be unity. The only unknown parameter in these ratios formed is the permittivity of the snow. The surface radar response is directly modeled using an empirical expression which is based on predictions of the IEM surface scattering model [37]computed over the range of surface parameters expected for snow. A correlation term is formed as for the volume case; once again, the correlation factor is considered to be unity. The unknowns in the empirical surface scattering models are the permittivity of the snow and a 18

general surface roughness term. Two equations are formed involving the surface scattering responses and the ratios of the volume responses. These equations are combined into a single equation, Eq. (22) in [8] which is only a function of the snow permittivity, Fs- The inversion algorithm amounts to finding the value of Ce which most nearly satisfies this equation. 4.3.1 Application to Measured Diurnal Data We first attempted to treat the Cadillac data set with this algorithm. The physical description of the snowpack in the Cadillac experiment seemed to agree generally with the inherent assumptions of the model, namely m, effects manifesting themselves primarily in the upper regions of the snowpack. Only the C-band data was used in the algorithm evaluation, which consists of testing a range of values of Cs to find the value most nearly satisfying Eq. (22) in [8]. We used a range which corresponded to values of this parameter which might reasonably be associated with snow: 1.1 < es < 5.0. The values of es which were the outputs of the algorithm with respect to the Cadillac data are shown in Figure 15. As seen, most of the twenty-nine separate sets of polarimetric radar measurements resulted in estimates of es which are within the bounds of "reasonable" results which we set. Only one data set produced a value pegged at the top of the range (eC = 5.0) and six were pegged at the bottom (e, = 5.0). Also shown in the figure is the nominal value of the permittivity of the top layer of the snowpack, based on ground truth measurements. The average measured density of the top layer was found to be about 0.5 g cm-3; such a density would correspond to a dry snow dielectric constant of 2.03, which is the nominal value shown in the figure. As described in [18], the permittivity of wet snow is modeled as a dry snow value given by, s = 1 + 1.7pds + 0.7ds, (13) which depends on density, Pds alone plus an incremental increase, which depends on snow wetness (my) alone. This incremental increase in the permittivity is given by [38], AC' -0.02m015 + 0.073ml.31 Ae, = 0.02m 1 + (/f)2 (14) 19

where f, is the frequency at which the permittivity is measured and fw =9.07 GHz is the relaxation frequency of water at 00C, and m, is expressed in percent. If the incremental increase is known, an accurate estimate of m, may be obtained by inverting Eq. (14). The incremental increase itself may really only be known if the density of the snow is known also, so that the base contribution to es represented by the dry snow component eds may be accounted for. Then the incremental increase is given by, Acs(mF ) = cs(mvp -cds(p) (15) This is one shortcoming of the algorithm by [8], and one which is not addressed by the authors, namely, the relatively large errors which result in trying to estimate mv of snow from a measurement of es alone, with no knowledge of the snow density [30]. For the purpose of evaluating the present algorithm however, we will compute estimates of mv for the Cadillac data set by subtracting the value of 2.03 (based on the measured density of 0.5 g cm-3) from each of the algorithm-estimated permittivities shown in Figure 15. The remainder of this operation is the incremental increase in the permitivity due to water, and this quantity can be inverted to obtain an estimate of m,. These estimates are shown in Figure 16, along with the values of my which were measured in the upper most layer as a function of time using the Snow Probe. All cases where the algorithm estimated value of es are less than the nominal dry-snow value of 2.03 are considered as having m, = 0. It is seen that the algorithm gives a reasonable performance in terms of its predictions that there were very high snow wetness levels present. About half of the cases examined result in estimates of m, = 0. Of the cases which give non-zero estimates of ma, very high wetness levels are indicated-similar to but in general exceeding the actual measurementswhich may be seen to-very roughly-follow the trends, as a function of time, which were observed in the measurements. As an additional test, we applied the algorithm to the Brighton diurnal data set, for which, as was seen in Figure 18, the top layer of the snowpack had relatively low levels of liquid water content. The values of es which we obtained from the application of the inversion algorithm to this data set are shown in Figure 17. The "nominal value" of Es indicated in this figure is based on the measured density of about 0.25 g cm-3 in the uppermost layer of the snowpack. As can be seen, most of the estimated values are lower than this nominal 20

level, but not by much. In particular, none of the nine separate cases treated resulted in permittivity estimates outside of the range (1.1 Es < 5.0) which are "reasonable". Also, the relatively very low estimates relative to those shown for the Cadillac case in Figure 15 demonstrate the algorithm is genuinely sensitive to this parameter. The associated estimates of m, generated in the same manner as was described for the Cadillac case above, that is, using our knowledge of the actual snow density in the topmost layer of the snowpack, are shown in Figure 18. Only two of the nine cases examined give non-zero estimates of the snow liquid water content. The actual measured values however are relatively low, with the possible exception of the values measured at about 6 p.m. (18.00 hours) which are seen to be above 5% liquid water content. This condition itself was a rather anomalous one, occurring as it did just prior due the refreezing of the snowpack top surface. 4.3.2 General Comments on the Wetness Retrieval Algorithm The central concept behind this algorithm appears to be that the relationships among the quantities oa,, (hh, and 0vvhh (the latter of which is related to the correlation between the co-polarized responses) for volume scattering can be distinguished from the corresponding relationships of these quantities associated with surface scattering. That this should be the case is not obvious. With respect to the co-polarized responses, for both surface and volume scattering it is generally the case that the vv-polarized response exceeds the hh- response increasingly more as incidence angle increases. With respect to the quantity avvhh, it has previously been noted by both [7] (the authors of the inversion algorithm) and others [4] 1 that the correlation coefficient for rough surface scattering is approximately unity. In addition, in the wetness inversion algorithm, the authors implicitly use a correlation coefficient of unity for the volume scattering as well. Thus it is difficult to see how the polarization relationships between the polarimetric quantities in volume and surface scattering respectively could allow discrimination between these two scattering mechanisms. Despite these questions about the fundamental concepts upon which the algorithm is constructed, it must be said that its performance with respect to the two quite different data sets of this present study - the Cadillac and Brighton diurnal experiments - is fairly impressive. Using only polarimetric radar data, the algorithm gave estimates of dielectric constant which were fairly comparable to those directly measured. In addi1 Specifically addresses the coefficient of variation for rough surfaces which, as we have indicated, carries essentially identical information as the correlation coefficient. 21

tion, when provided with density values for the two respective snowpacks, the algorithm was able to produce estimates of liquid water content reasonably close to those found to be present by direct measurements. 4.3.3 On the Question of the Behavior of o~ with Increasing m, The authors of this wetness inversion algorithm make a case for the increase in the backscatter levels for wet snow. While it is true that such a phenomenon is readily predicted theoretically, it has only been observed experimentally in some exceptional cases. The Cadillac scenario might be considered such a case, where a combination of circumstances (ice lens near the surface and exceptionally warm weather) produced, to a small degree, a surface scattering effect. In general however, the depression of the radar backscatter level has been almost universally observed in practice [12, 13, 10, 15, 23] at both microwave (including the present study) and millimeter wave frequencies. In [12] examples are given of the opposite case, of increasing ~, but these cases were brought on by the presence of rain in a time frame very close to the measurements in one case, the presence of hail in another case, and in one other case, wet snow which was roughened with a shovel. The rain scenario, in fact, mirrors the circumstances under which the experiment in the Oztal Alps was conducted-an experiment which produced the data in which the inversion algorithm was originally based. The preponderance of evidence with respect to a decreased backscatter suggests that models of wet snow may not be realistic. Possible explanations for the departure from theoretical predictions may be that (i) liquid water drains from the upper levels of the snowpack before it becomes abundant enough to increase the dielectric constant significantly or, as may well be the case with the Cadillac data (ii) drainage from the immediate top of the layer helps to create a sort of matching layer leading to small surface scattering even for high wetness levels near the top surface. 5 Summary This paper has described the results of two types of polarimetric radar experiments which were carried out on snowpacks. In Section 3 we presented results and analysis for measurements which were made at C-, and X-band on the bare ground and then three successively deeper (20,60, and 102 cm) layers of artificial dry snow. The 22

details of the physical character of the snowpack and the environmental conditions associated with the experiments made these results especially amenable to comparison with discrete -particle -based theoretical modeling techniques. It was shown, however, that these techniques did not give reasonable agreement with the experimental observations. A subsequent analysis of the data was presented by which polarimetric scattering and extinction quantities intrinsic to the snow medium were retrieved. In Section 4, results for backscatter collected during partial diurnal cycles were presented, along with the results of extensive measurements of snow liquid water content which were made concurrently. This data was used in an attempt to confirm the validity of an algorithm which has recently been developed by [8] for the retrieval of snow liquid water content from polarimetric C-band measurements. Although an examination of the conceptual framework of the algorithm reveals certain basic assumptions which seem difficult to justify, its performance with respect to the two separate diurnal data sets is fairly impressive. The algorithm was able to correctly characterize the Cadillac and Brighton snowpacks (top layer) as very wet and reasonably dry, respectively. 23

References [1] J. Koskinen, L. Kurvonen, V. Jaaskelainen, and M. Hallikainen, "Capability of radar and microwave radiometer to classify snow types in forested areas," in Proc. IGARSS '94, p. 1283, 1994. Vol. II. [2] F. T. Ulaby, P. Siqueira, and K. Sarabandi, "A hybrid electromagnetic-statistical approach for characterizing mmw scattering by terrain," in AGARD Conf. Proc. 501, 1993. [3] R. M. Narayanan and S. R. Jackson, "Snow cover classification using millimeter-wave radar imagery," in Proc. IGARSS '94, 1994. Vol. IV. [4] J. J. van Zyl, H. A. Zebker, and C. Elachi, "Imaging radar polarization signatures: Theory and observation," Radio Science, vol. 22, no. 4, pp. 529-543, 1987. [5] J. J. van Zyl, "Unsupervised classification of scattering behavior using radar polarimetry data," IEEE Trans. Geosci. Rem. Sens., vol. 27, no. 1, pp. 36-45, 1989. [6] D. L. Evans, T. G. Farr, J. J. van Zyl, and H. A. Zebker, "Radar polarimetry: Analysis tools and applications," IEEE Trans. Geosci. Rem. Sens., vol. 26, no. 6, pp. 774-789, 1988. [7] J. Shi and J. Dozier, "Radar backscattering response to wet snow," in Proc. IGARSS '92, p. 927, 1992. Vol. II. [8] J. C. Shi and J. Dozier, "Inferring snow wetness using c-band data from sir-c's polarimetric synthetic aperture radar," IEEE Trans. Geosci. Rem. Sens., vol. 33, no. 4, pp. 905-914, 1995. [9] D. T. Davis, Z. Chen, L. Tsang, J. N. Hwang, and A. T. C. Chang, "Retrieval of snow parameters by iterative inversion of a neural network," IEEE Trans. Geosci. Rem. Sens., vol. 31, no. 4, pp. 842-852, 1993. [10] W. H. Stiles and F T..Ulaby, The active and passive microwave response to snow parameters 1. wetness," J. Geophys. Res., vol. 85, no. C2, pp. 1037-1044, 1980. [11] F. T. Ulaby and W. H. Stiles, "The active and passive microwave response to snow parameters 2. water equivalent of dry snow," J. Geophys. Res., vol. 85, no. C2, pp. 1045-1049, 1980. 24

[12] C. Matzler, E. Schanda, and W. Good, 'Towards the definition of optimum sensor specifications for microwave remote sensing of snow," IEEE Trans. Geosci. Rem. Sens., vol. GE-20, no. 1, 1982. [13] C. Matzler and E. Schanda, "Snow mapping with active microwave sensors," Int. J. Rem. Sens., vol. 5, no. 2, pp.409-422,1984. [14] A. T. C. Chang, J. L. Foster, M. Owe, and D. K. Hall, "Passive and active microwave studies of wet snowpack properties," Nordic Hydrology, vol. 16, pp. 57-66, 1985. [15] F. T. Ulaby, T. F. Haddock, R. T. Austin, and Y. Kuga, "Millimeter-wave radar scattering from snow: 2. comparison of theory with experimental observations," Radio Science, vol. 26, no. 2, pp. 343-351, 1991. [16] J. B. Mead, P. S. Chang, S. P. Lohmeier, P. M. Langlois, and R. McIntosh, "Polarimetric observations and theory of millimeter-wave backscatter from snow cover," IEEE Trans. Antennas Propagat., vol. 41, no. 1, pp. 38-46, 1993. [17] M. T. Hallikainen, F. T. Ulaby, and T. E. Van Deventer, "Extinction behavior of dry snow in the 18- to 90-ghz range," IEEE Trans. Geosci. Rem. Sens., vol. GE-25, no. 6, pp. 737-745, 1987. [18] J. R. Kendra, F. T. Ulaby, and K. Sarabandi, "Snow probe for in situ determination of wetness and density," IEEE Trans. Geosci. Rem. Sens., vol. 32, no. 6, pp. 1152-1159, 1994. [19] M. A. Tassoudji, K. Sarabandi, and F. T. Ulaby, "Design consideration and implementation of the lcx polarimetric scatterometer (polarscat)," Tech. Rep. Rep. 022486-T-2, The University of Michigan, Radiation Lab, 1989. [20] Y. Oh, K. Sarabandi, and F. T. Ulaby, "An empirical model and an inversion technique for radar scattering from bare soil surfaces," IEEE Trans. Geosci. Rem. Sens., vol. 30, no. 2, pp. 370-381, 1992. [21] K. Sarabandi, Y Oh, and T., Ulaby, "Measurement and calibration of differential mueller matrix of distributed targets," IEEE Trans. Antennas Propagat., vol. 40, pp. 1524-1532, Dec. 1992. [22] D. R. Brunfeldt, '"Theory and design of a field-portable dielectric measurement system," in Proc. IGARSS '87, pp. 559-563,1987. Vol. I. 25

[23] P. Chang, J. Mead, R. McIntosh, R. Davis, and H. Boyne, "A detailed study of the backscatter characteristics of snowcover measured at 35, 95, and 225 ghz," in Proc. IGARSS '94, p. 1932, 1994. Vol. IV. [24] M. T. Hallikainen, F. T. Ulaby, M. C. Dobson, M. A. El-Rayes, and L. Wu, "Microwave dielectric behavior of wet soil -part i: Empirical models and experimental observations," IEEE Trans. Geosci. Rem. Sens., vol. GE-23, pp. 25-34, 1985. [25] A. Ishimaru and Y. Kuga, "Attenuation constant of a coherent field in a dense distribution of particles," J. Opt. Soc. Am., vol. 72, pp. 1317-1320, 1982. [26] L. Tsang and J. A. Kong, "Scattering of electromagnetic waves from a dense medium consisting of correlated mie scatterers with size distributions and applications to dry snow," Journal of Electromagnetic Waves and Applications, vol. 6, no. 3, pp. 265-296, 1992. [27] E. Nyfors, "On the dielectric properties of dry snow in the 800 mhz to 13 ghz range," Tech. Rep. S135, Helsinki University of Technology, Radio Laboratory, 1982. [28] Y. Kuga, F. T. Ulaby, T. F. Haddock, and R. D. DeRoo, "Millimeter-wave radar scattering from snow: 1. radiative transfer model," Radio Science, vol. 26, no. 2, pp. 329-341, 1991. [29] L. M. Zurk, K. H. Ding, L. Tsang, and D. P. Winebrenner, "Monte carlo simulations of the extinction rate of densely packed spheres with clustered and non-clustered geometries," in Proc. IGARSS '94, pp. Vol. I, p. 535, 1994. [30] J. R. Kendra, Microwave Remote Sensing of Snow: An Empirical I Theoretical Scattering Modelfor Dense Random Media. PhD thesis, The University of Michigan, 1995. [31] F. T. Ulaby and C. Elachi, eds., Radar Polarimetryfor Geoscience Applications. Dedham, MA: Artech House, Inc., 1990. [32] E. P. W. Attema,, and F. T. Ulaby, "Vegetation modeled as a water cloud," Radio Science, vol. 13, pp. 357-364, 1978. [33] Y. Oh, K. Sarabandi, and F. T. Ulaby, "An empirical model for phase difference statistics of rough surfaces," in Proc. IGARSS '93, 1993. 26

[34] K. Sarabandi, "Derivation of phase statistics from the mueller matrix," Radio Science, vol. 27, no. 5, pp. 553-560, 1992. [35] J. L. Zhou and A. L. Tits, "User's guide for fsqp version 3.3b: A fortran code for solving constrained nonlinear (minimax) optimization problems, generating iterates satisfying all inequality and linear constraints," Tech. Rep. TR-92-107r3, Systems Research Center, The University of Maryland, September 1993. [36] L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing. New York, NY: Wiley Interscience, 1985. [37] A. K. Fung, Z. Li, and K. S. Chen, "Backscatttering from a randomly rough dielectric surface," IEEE Trans. Geosci. Rem. Sens., vol. 30, no. 2, pp. 356-369,1992. [38] M. Hallikainen, F. T. Ulaby, and M. Abdelrazik, "Dielectric properties of snow in the 3 to 37 ghz range," IEEE Trans. Antennas Propagat., vol. AP-34, pp. 1329-1339, 1986. 27

List of Figures 1 Temperature variations during the course of the radar experiments on artificial snow..... 30 2 Photograph of artificial snow particles. Major divisions of ruler shown are millimeters.... 31 3 Measured (discrete marks) and predicted (lines) backscatter from bare ground......... 32 4 Backscatter measurements (VV and VH) made of three different snow (artificial) depths at C-band and X-band. Also shown for reference are simulations of the pure ground scattering expected and also the contributions from the ground and the snow top surface......... 33 5 Comparison of measured backscatter results for dry (artifical) snow at C- and X-band with optimal RT predictions obtained by treating the particle size as a free parameter. Optimal snow particle diameters are (1) C-band: 1.7 mm, and (2) X-band: 0.9 mm. Measured snow particle diameter is 0.27 ~ 0.11 mm............................... 34 6 First-order volume scattering mechanisms in a layer of scatterers................ 35 7 Comparison of snow data to radiative transfer-type model computations. Parameters for model are obtained through an optimization process which uses this data............... 36 8 Comparison of C- and X-band snow data to radiative transfer-type model predictions. Model is parameterized though analysis of 400 data.......................... 37 9 Comparison of C- and X-band snow data to radiative transfer-type model predictions, for a 102-cm layer. Model is parameterized though analysis of 40~ data............... 38 10 Phase statistics for the 60 cm artifical snow layer. Shown are the measured values and those calculated using an RT model with associated parameters given in Table 2.......... 38 11 Co- and cross-polarized backscatter results from Brighton (partial) diurnal experiment. Incidence angle is 400....................................... 39 12 Vertical profiles of snow wetness measured as a function of time. Measurements were taken with the snow probe during the Brighton diurnal experiment.................. 39 13 Co- and cross-polarized backscatter results from Cadillac (partial) diurnal experiment. Incidence angle is 40....................................................... 40 14 Vertical profiles of snow wetness measured as a function of time. Measurements were taken with the snow probe during the Cadillac diurnal experiment.................. 40 28

15 Application of inversion algorithm for snow wetness applied to Cadillac diurnal data set. Shown is the actual output from the algorithm, snow permittivity E. The "nominal level" shown is based on the expected dry snow permittivity based on ground truth measurements of density............................................. 41 16 Snow wetness inversion algorithm results from Cadillac diurnal data set compared with actual measured values of mv in the uppermost layer of the snowpack.............. 41 17 Application of inversion algorithm for snow wetness applied to Brighton diurnal data set. Shown is the actual output from the algorithm, snow permittivity ae. The "nominal level" shown is based on the expected dry snow permittivity based on ground truth measurements of density............................................. 42 18 Snow wetness inversion algorithm results from Brighton diurnal data set compared with actual measured values of m, in the uppermost layer of the snowpack.............. 42 29

TEMPERATURE DURING BRIGHTON EXPMT 10.0 5.0 0 -10.0 H H LO 60 CM -15.0 -20.0 BARE 102 CM 20 CM.25.0 I — 16 21 26 3 8 (FEB) (MAR) DAY Figure 1: Temperature variations during the course of the radar experiments on artificial snow 30

I I I I: Figure 2: Photograph of artificial snow particles. Major divisions of ruler shown are millimeters. 31

0.0 -5.0 -10.0 -15.0 -20.0 % -25.0 -30.0 -35.0 -40.0 -45.0 CAy 11.... * I. ' ' '.' I ' ' '.' I.... — E — HH ------ VH 6 ---- e — ------ - - " - - - 0 A A ~' ---^.....,I -:11I.. 10.0 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (a) C-band l l..,........................... I.. -5.0 -10.0 -15.0 -20.0 ~t -25.0 -30.0 -35.0 -40.0 -45.0 --—. vv — o — VV - -- - HH ------ VH.- - - - I - - -....*I -;- * I... - ~ A -'.^ i -. I i -CA IlA-1,I.... I.... I.... I........ - 10.V 10.0 20.0 30.0 40.0 50.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (b) X-band Figure 3: Measured (discrete marks) and predicted (lines) backscatter from bare ground. 32

-5.0 I I I I I I I I -15.0 F - ID.. ----- ~j- — 0_ "O -25.0 - 8 8 8 VH 102 cm 60 cm ---- 20 cm 0 gnd + top, sim. 0 gnd, sim. -35.0 - -45.0 10.0 8 8 I, I I I I 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (deg.) (a) C-band -5.0 -15.0 1-s - o 0 0 -25.0 - " ^ " - | - vv" - ~ — I 0 VH E —3 102cm 60 cm ------ 20 cm o gnd + top, sim. E gnd, sim. -35.0 -45.0 10. 0 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (deg.) (b) X-band Figure 4: Backscatter measurements (VV and VH) made of three different snow (artificial) depths at C-band and X-band. Also shown for reference are simulations of the pure ground scattering expected and also the contributions from the ground and the snow top surface. 33

.20.0 co 0 -30.0 (a) C-band measured. ------ -.16.0 '" 200 30.0 0. -360 10.0. 1.0 incidence prosic (deg-3 imulatio. (b) C-band optimal.T \ ------------- w A -16.0 \ \ 60 cm \ I DO nl I 20.0, I ncide e, e (deg.) )X-band optimal T simulatio (c) X band measured. -t C and ith optmale s C-band-ad. '.9 mm 34

Figure 6: First-order volume scattering mechanisms in a layer of scatterers. 35

0.0 -10.0 0 ~& -20.0 -30.0 -40.0 0.0 20.0 40.0 60.0 80.0 100.0 Depth (cm) (a) C-band, 400 0.0 -10.0 "a - -20.0 -30.0 -40.0 0.0 20.0 40.0 60.0 80.0 100.0 Depth (cm) (b) X-band, 40~ Figure 7: Comparison of snow data to radiative transfer-type model computations. Parameters for model are obtained through an optimization process which uses this data. 36

0.0 -10.0 0 -o 0 -20.0 —,~ — C-VV — C — X-VH.-~ —. C-VH o o ---.... ---" " " "A" — ---- n --- —. a G, -—................................... -30.0 - -40.0 10. 0 0O............................ 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (deg.) (a) 20-cm layer 0.0 -10.0 -20.0 -30.0.. I. I.. ~ I I I I. I. I I.., 1 I ~a 0 10 — o-<~.X-VV —,- - C-VV — < ---- X-VH - ---- -........-...... C-VH -----------------— A. *.......... A............................. 20.0 30.0 40.0 50.0 60.0 70.0 20.0 30.0 40.0 50.0 60.0 70.0 -40.0 L10.0 Incidence Angle (deg.) (b) 60-cm layer Figure 8: Comparison of C- and X-band snow data to radiative transfer-type model predictions. Model is parameterized though analysis of 400 data. 37

0.0 -10.0 I............................ 0 tD -20.0 I.....* 1..... 1 -..r.............. 1 ' ' ' ----—. X-VV - 49 - - C-VV ----- X-VH ------------- --------- C-VH — o < —o ----— D-.& -_ o................. I...............I L.......I...J................................................... -30.0 F Af A A -4U.U. 10.0 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (deg.) Figure 9: Comparison of C- and X-band snow data to radiative transfer-type model predictions, for a 102-cm layer. Model is parameterized though analysis of 40~ data. 0 o t_ U o '3 1.00 0.80 0.60 0.40 0.20 0 0 0 D'-" --- —B- - - -. - -.. '-......a.."'' '-.................. 0 0 MEAS. --—.. --- HYB. RT. MODEL h C-/ C3 0 v.ua o.N 200.0 190.0 180.0 170.0 160.0................-. I-.-..-. B ---. -C.C-u-.C 0 0 MEAS. --— 3 ---- HYB. RT. MODEL I O.o I....................... I........ I.... I.... I... I 10.0 20.0 30.0 40.0 50.0 Incidence Angle (deg.) (a) Degree of Correlation 60.0 70.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (deg.) (b) Co-polarized Phase Difference Figure 10: Phase statistics for the 60 cm artifical snow layer. Shown are the measured values and those calculated using an RT model with associated parameters given in Table 2. 38

0.0 -10.0 -20.0 -30.0 -— ~.... h — -— e --- hh -40.0 ---- vh — v — hv -50.0........ I 10.0 12.0 14.0 16.0 18.0 20.0 Time of Day (a) C-band 0.0 O ~ ~,......-..,,. -10.0 -20.0 -30.0 40.0 a., A -— G —. hh ----- vh -- " — hv -50.0 1... I 10.0 12.0.. I... I........ _ 14.0 16.0 Time of Day (b) X-band 18.0 20.0 Figure 11: Co- and cross-polarized backscatter results from Brighton (partial) diurnal experiment. Incidence angle is 40~. ' L 15 0-e to 0 3.-:O ' O I..-:.....:..:..:... A.: -... i-:.:..: ---. An. ' I 020 20 20 18 40 16 60D 14 80 12 100 10 Time of Day (hr) Height above ground (cm) Figure 12: Vertical profiles of snow wetness measured as a function of time. Measurements were taken with the snow probe during the Brighton diurnal experiment. 39

0..1 -10.0 --— O- VV.... --- HHH -— A --- VH -- --- HV — 4 -~-v Q ' $"^c '^"^WB —<^^p~ -20.0 -30.0 k 0.. ' I.. -- - I ' ' ' I ' ' ' I — 0 --- VV -— Q — HH -10.0 -— & --- VH -- - HV -20.0 -30.0 - ': 8 ~ — -40.0 -50.0............. 9.50 10.50 11.50 12.50 13.50 14.50 15.50 16.50 Time (hrs) (b) X-band -40.0 -50.0....... 9.50 10.50 11.50 12.50 13.50 14.50 15.50 16.50 Time (hrs) (a) C-band Figure 13: Co- and cross-polarized backscatter results from Cadillac (partial) diurnal experiment. Incidence angle is 40~.. A *a '' ''' ` I 1 '' 16 0 15 5 14 10 13 15 12 20 11 25 10 Time of Day (hr) Height above ground (cm) Figure 14: Vertical profiles of snow wetness measured as a function of time. Measurements were taken with the snow probe during the Cadillac diurnal experiment. 40

8.0 7.0 6.0 5.0 CA 0 co 40 0 r) Co) 4.0 ~ I I ' I ~ ~ UPPER LIMIT (e's=5.0).-. --- NOMINAL LEVEL (~'s=2.03) -LOWER LIMIT (E's= 1.1) - - - - SHI-DOZIER ALGORITHM (1995), - *. *;', ' * ' ' ' ' ' '. ', I ' [ ---- --- -- -* --- - —; - - -— / —y- - -* —7 T-E -- i- -- *' ---- - ----- - - -- --- - — T --- do - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -I-.I.-I-. I. --- —--— _ 3.0 2.0 1.0 An V.-X 9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 Time of Day (hr) Figure 15: Application of inversion algorithm for snow wetness applied to Cadillac diurnal data set. Shown is the actual output from the algorithm, snow permittivity er. The "nominal level" shown is based on the expected dry snow permittivity based on ground truth measurements of density. c-t C3 I CO rA 19 17 15 13 11 9 7 5 3 1 -1 9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 Time of Day (hr) Figure 16: Snow wetness inversion algorithm results from Cadillac diurnal data set compared with actual measured values of m, in the uppermost layer of the snowpack. 41

3.0 0) o 0 i U o o c) 2.0 I ' I I ' I * I UPPER LIMIT (~='5.0) NOMINAL LEVEL (E'= 1.47) LOWER LIMIT (E's=1.1) - -. - - SHI-DOZIER ALGORITHM (1995),e -- ------ ---.' --------------------- -------------------------------- I. I. I. I. I, I 1.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0 22.0 Time of Day (hr) Figure 17: Application of inversion algorithm for snow wetness applied to Brighton diurnal data set. Shown is the actual output from the algorithm, snow permittivity ea. The "nominal level" shown is based on the expected dry snow permittivity based on ground truth measurements of density. E I co u: C3 <) On 0 cl) 19 17 15 13 11 9 7 5 3 1 -1 8 10 12 14 16 18 20 22 Time of Day (hr) Figure 18: Snow wetness inversion algorithm results from Brighton diurnal data set compared with actual measured values of m, in the uppermost layer of the snowpack. 42

Appendix B A Hybrid Experimental/Theoretical Model for Dense Random Media

A Hybrid Experimental / Theoretical Scattering Model for Dense Random Media John R. Kendra and Kamal Sarabandi Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Abstract The subject of scattering of electromagnetic waves by dense media has been one of intense interest in recent years. The present paper describes polarimetric backscatter measurements made at Ku-band on layers of a dense medium under very carefully controlled circumstances. The experiments have a dual purpose: (i) to evaluated the degree to which the experimental observations are predicted by theoretical, particle-based, random media models; and (ii) to test a proposed hybrid model by which the scattering and extinction properties of a dense medium are characterized experimentally, allowing future modeling of the polarimetric response for any arbitrary configuration of the medium. The hybrid model assumes that first-order vector radiative transfer is a suitable theoretical structure providing the extinction and phase matrix components are appropriately specified; the specification is accomplished through an inversion algorithm involving polarimetric backscatter measurements. The major conclusions of the study are that: (i) the hybrid model is an adequate description of the dense medium scattering behavior; (ii) conventional radiative transfer (RT) appears to give a reasonable estimate of the observed radar response, but dense medium RT gives a very low estimate; and finally, (iii) the phase function of the effective volume scattering element of the medium, obtained via the hybrid model, suggests a larger effective scatterer than the physical ones. 1

1 Introduction Evidence has recently been presented [1, 2] which suggests that existing particle-based theories are inadequate for modeling very dense random media. This raises the very fundamental question of what recourse exists for simulating the effects of a dense medium, for example, for the case of a snow layer over sea ice, or on a forest floor. The ability to generate realistic predictions is an essential prerequisite to any scheme for retrieving physical characteristics about a target from remotely sensed data. If the behavior of densely packed discrete particles is difficult to model, it is evident that the issue is only exacerbated when dealing with materials of a very complex physical character, for example the amorphous interconnected form which snowpacks can often-or even usually-take. There are of course theoretical techniques which could be applied to such cases; these are field-based techniques in which the medium is described as a fluctuating dielectric constant. To obtain a solution, the Born approximation or the Distorted Born approximation is usually applied. In practice however, although some studies have investigated this technique [3] the particular material characteristic which is required as an input, the correlation function of the medium, is exceedingly difficult to obtain. The standard technique for its measurement is the "thin sections" technique. In this technique [3, 4] a super-cooled liquid is allowed to fill the pore spaces of a snow sample. After freezing the samples, they are shaved on a microtome and polished, treated with a contrast enhancer, and digitized. Not only is this a very arduous process, but it has been recently shown [5] that the correlation function must be known with high accuracy, including its tail region, to obtain accurate prediction of scattering; that is, it is not sufficient to know that the function approximates, say, a gaussian or exponential function. In additon to these obstacles, it is also true that the distorted Born approximation does not account for multiple scattering of the incoherent wave [6]. Apart from the theoretical approaches addressed above, purely empirical approaches may be considered; these however have the obvious limitation that the entire parameter space of the target cannot be sufficiently well known to allow estimation of more specific target properties. To circumvent the difficulties associated with the above mentioned techniques and to offer some means by which realistic modeling of dense media might be accomplished, a new hybrid experimental / theoretical modeling scheme is introduced in this paper. In the sections that follow, we will describe the hybrid model, and address its validity with actual polarimetric backscatter measruements on a dense medium. The measured 2

results will also be compared with certain discrete-particle based random media theories. 2 Hybrid Model Concept The hybrid model involves two major assumptions. The first is that radiative transfer (RT) theory is applicable. That is, that the flow of electromagnetic energy through any dense medium obeys the fundamental equation of transfer, dI(~r) -icc(s)I( s r) + j f T(3, v)I(st, r)di' (1) ds -0 =0 where I is the 4 x 1 Stokes vector, rCC is the 4 x 4 extinction matrix, and P? is the 4 x 4 phase matrix. The phase and extinction matrices for a medium may be formulated in terms of the characteristics of a single particle if the medium consists of uniform particles or a weighted average of the individual characteristics of various particles if the medium is heterogeneous [7]. This type of formulation, known as conventional radiative transfer (CRT), contains two simplifying assumptions: (i) that all particles are in the far-field of one another and (ii) that the particle positions are completely uncorrelated with one another, thereby eliminating any coherence effects. It has also been shown [8] that for Rayleigh particles, a rigorous field-based approach to scattering in a dense medium, where particle correlations are considerd and exact wave transformations from one particle to another are employed, yields a solution which may be couched in exactly the same form as (1). It was found from this analysis that the phase function is the Rayleigh phase function, similar to CRT. However, the scale of the phase function, which is provided by the scattering albedo, and the extinction differ from CRT. Thus the assumption that such an equation describes the mechanism is reasonable, although it has not been shown to hold for the case of a rigorous field-based analysis of non-Rayleigh particles [9]. We have already mentioned above studies in which evidence is provided suggesting that the available techniques by which Cc and 1P are computed do not appear to be adequate. In addition, for complex or amorphous materials, where a "particle" cannot be unambiguously identified, it is not possible to formulate these quantities. Therefore, in the present hybrid model, no attempt is made to specify them based on the physical characteristics of the medium. Instead, a direct measurement process will be carried out, effectively interrogating the medium, and from these measurements an inversion process will be used to retrieve ice and P. 3

In order to make this approach tractable, a second major assumption is required. This is that, not only does radiative transfer apply, but that the first-order evaluation of it is sufficient. The first-order solution to the radiative transfer equation for the case of an isotropic layer with smooth surfaces above and below has been derived in [7]. For, IS(,o ) = Lm(Po)II(-o ) (2) where Lm is transformation matrix related to the Mueller matrix Mf [7] by, 1 Lm Acos (3) A cos Oj where A is the illuminated area, and ei the angle of incidence, we have, Lm(1o) = - 21 [Y.23PbsR23 + d T23Pbi + -Pbi23 + Pbs] T12 (4) o L L where, 1 - exp(-2rced/go) Y - 2Le/) (5) L = exp(2rced/ )g, (6) where ie is the extinction (a scalar quantity) and where Pbs and Tbi are the backscatter and bistatic components of the phase matrix, respectively; HRpq and Tpq are the respective Fresnel reflectivity and transmissivity matrices associated with intensity propagatg in m m g tiar mediumm m q, and. is the cosine of the refracted angle in the medium. The four scattering mechanisms described by Eq. (4) are shown in Figure 1. The scattering event depicted in terms (A) and (B) are associated with the backscatter component of the phase function,!bs; terms (C) and (D) are associated with the bistatic component of the phase function, Tbi. In the figure, we have intentionally represented the scattering elements as clusters to underscore the point that, in this treatment, we are considering effective particles, which may comprise correlated groups of individual physical particles 4

and/or multiple scattering effects. 3 Rayleigh Model for Phase and Extinction Matrices In this section, we examine, as an aid to inversion of the hybrid model, the essential character of the extinction matrix Kc (analysis begins with the most general case of a matrix extinction as derived in [10]), and the two components of the phase matrix Pbs and Pbi, for the special case of scattering by a layer of Rayleigh particles. This scenario corresponds to the isotropic case mentioned above, where IC and Tbs are constant matrices, but Pbi is a function of angle. The analysis is motivated by the intuition that the number of measurements required and the complexity of the ensuing inversion operation could be reduced through knowledge of the general form of the unknown matrices; that is, knowledge of which elements are non-zero, which are independent, and the nature of the dependences between elements. As a starting point, we consider the effective particles of the layer (as depicted in Figure 1) to be Rayleigh particles. As we will show in the following sections, this will allow us to derive symbolically the specific mathematical structure of the quantities of interest - the extinction matrix, and certain components of the phase matrix consistant with a first-order approach - in terms of the elements of the polarizability tensor of the particles. This derivation follows closely that found in [11], except that in that case only axisymmetric particles were considered, whereas, in the present case, no assumptions are made about the particle's symmetry properties. 3.1 Rayleigh Theory: Scattering Matrix Element Representation For particles with dimensions small relative to the wavelength, scattering may be described in terms of a polarizability tensor. A scattered electric field vector is given by [12]: ES= -0 eikrs Xis xp (7) 4nr where fk is the unit vector in the direction of propagation of the scattered wave, and p is the induced dipole moment determined from the solution of the Laplace equation. If the induced dipole moments (px, Py, Pz) 5

are derived for (x, y, z) polarized incident fields, then The scattered field may be expressed as ES = -)eikr s x ks x (D. Ei) (8) 47Cr where, D is the 3 x 3 polarizability tensor having as its columns the vectors Px, Py, and pz, and E1 is the incident field vector, [Ex, Ey,Ez]T. The polarizability tensor D is a function of the geometry and dielectric constant of the Rayleigh particle. In general, D will be a function of the orientation of the particle. Since, eik~r Es = eir Ei (9) where, S (V vhI (10) k\Shv Shh/ it can be shown [13] that the scattering matrix elements corresponding to the scattered polarization s and the incident polarization qi such that Ps E {sfis} and qi E {vi,fii}, are given by Spq = 4 Ps. (D. i).() 3.2 Derivation of Extinction and Phase Matrix Elements The extinction and phase matrix elements are expressed in terms of the ensemble-averaged quantities (Spq) [10] and (SpqS* p) [7], respectively. In this section, we present the derivation of the forms of these elements in terms of the elements of the polarizability tensor in global coordinates. 3.2.1 Extinction Matrix Elements We consider first terms of the form (Spq). From (11), (Spq) = ' (s - (ID- qi)). (12) 47t P 6

Since the polarization vectors themselves are not a function of the particle orientation, the averaging process is performed only on!D itself: k2 (Spq) = 4 s(()q). (13) The form of!D in local coordinates is transformed from that in global coordinates (!l) through the application of transformation matrices involving the Eulerian rotation angles a, P, and y(see, for example, [14], pp. 158 -160): DA AT. V -A (14) where the elements of the 3 x 3 transformation matrix A are trigonometric functions of a, 3, and y, e.g., all = cosYcoscosa- sinysina a12 = cosycos psina + sinycosa (15) al2 = -cosysinof Therefore for each element of D, the polarizability tensor in local coordinates is: 3 di (JAT. V.A)ij = - (amianj)dmn. (16) m,n= 1 It is seen from (16) that the symmetry property of ) is preserved across a coordinate transformation. For extinction, the quantity (Spq) is computed in the forward direction only. Because of the assumption of classical radiative transfer that the positions of particles are uncorrelated, (Spq) will be zero for every direction but forward scattering. Since we have assumed there is no orientation dependence, we can assign the simplest possible forms to the vertical and horizontal polarization vectors, incident and scattered, and compute (Spq) from (13). Considering an incident direction along the positive x-axis, using the Forward Scattering Convention (FSA) [7], 7

we have, hi = [0 1 0] = h (17) vi [OO-l] =T Ignoring the factor k, and applying (13), (Svv) = (d33) (18) (Shh) = (d2) (Shv) = -(d3) (Svh) = -(d2) = -(d3) From (16), 2ir X 21 3 (d1j) = / f/ I (amianj)dmnP(a)p(p) p(y)daddy (19) J=0 P=0o Jy=o m,n= where, for completely uniform orientation, the respective probability density functions are: 1 p(a)= 2 (20) (20) () sinf 2 P(Y) = 21 over the limits indicated on the integrals in (19). These integrals may be evaluated without difficulty and the results are: (d22) = (d + d2 + d33) (21) 321 (d33) = (d22) (d23) = 0. 8

Thus for the scenario under consideration, the extinction is found to be a scalar (as was already assumed in Eq. (4) and - when the result above is used in the definition [7] for the extinction elements - we obtain the following expression: Ke = -2Re[ 7n (dl+d22+d3)]' (22) 3ko 3 where no is the particle number-density. 3.2.2 Phase Matrix Elements The derivation of the phase matrix elements involves evaluation of quantities of the form (SpqS ). Referring back to Eq. (11), we have: k4 (SpqSmn) = (4;)2 ((Os D qj)(mis' i D nii)). (23) The derivation is similar to that performed for the extinction matrix elements but is considerably more involved. The problem of evaluating Eq. (23) essentially reduces to evaluating the spatial ensemble average of the Kronecker tensor product, {<D Z~ D), (24) where ZD is the polarizability tensor in local coordinates whose elements are given by Eq. (16). The Kronecker tensor product, however, inflates the formulation greatly since each element dij in D is composed of approximately twelve terms (since some of the elements in the transformation matrix A have two terms), each element in (Do 0 P*), a 9 x 9 matrix, has on the order of one hundred and forty-four terms, each requiring analytical integration with respect to the three Eulerian rotation angles a, [, and y. The details of this 9

derivation are given in [13]. The result of the analysis is the following extremely simple form: bs = P1 (P1 - P2) 0 0 O O 2(P1 -P2) P1 0 0 O O 0 0 -2(P +p2) 0 0 0 0 2 (PI - 3P2) (25) and, Pbi = P(l + cos2(20)) - P2sin2(2) (P1 -P2) 2(P1 -P2) PI 0 0 0 0 (26) 0 0 2 cos20(Pi +P2) 0 0 0 - cos 2(P - 3P2) where, 1 4 2 P1i= -A+ B+ C 5 15 15 P2= -A —B+-4C 15 15 15 A = I{d112 + 1412 + l 12 B 12+ 2+143 B = Rd1212 + Id'312 + d2312 ) Re(dll d2*2 + dld33 + d22d33 ). I \C"122TU1 33 Thus, an analytical evaluation of the theoretical problem of backscattering from a layer of general Rayleigh particles reveals that scattering and extinction within the layer is described by just three parameters: a scalar extinction, Ke and the two scalar parameters P1 and P2. This extremely simple form can potentially be exploited in the inversion process associated with the hybrid model concept. 10

This completes the theoretical framework of the Hybrid Model concept. The remaining elements are an experimental process and a numerical inversion algorithm. These are described in the following sections. 4 Experimental Process An experimental effort was carried out to test the validity of the Hybrid Model concept. A photograph of the experimental setup is shown in Figure 2. The experimental process was straightforward. Polarimetric backscatter measurements of a dense material were made with a Ku-band, network-analyzer-based radar. Details of the specific material examined are given in Section 6. The material examined was contained inside a large 1.8 x 1.8 meter wooden "sandbox". The sandbox itself was mounted on top of a turntable to allow for independent spatial measurements. The sides of the sandbox were two-tiered with the top tier removable so that it could have a maximum depth of either 23 cm or 43 cm. This was useful to prevent an excessive amount of "shadowing" of the surface of the layer for shallow depths. The radar was mounted on a manually adjustable, telescoping mechanical lift. The radar mount was such that manual adjustment of the incidence angle could be accomplished by simply tilting it and locking it into place. The range from the aperture of the radar antenna to the layer surface, along the boresite direction of the antenna was always 3 meters. This corresponds to approximately D2/X where D is the dimension of the square aperture of the antenna, and was sufficiently small to allow a full range of incidence angles to be observed in an indoor setting. As is described fully in [13], the radar consists of upconversion, downconversion, and polarization selection circuitry attached to a square horn antenna via an orthomode transducer. A coherent IF signal at 1-3 GHz is provided by the network analyzer and mixed with a local oscillator at the radar box to produce a 15 -17 GHz RF. The design and operation of the radar follows closely that of the LCX POLARSCAT system described in [15]. A compact geometry for the antenna is enabled through the use of a dielectric lens to correct the phase error at the aperture. The final beamwidth for the lens-corrected antenna (one-way, HPFW) was 5.3~ for both E- and H-plane. Two target parameters were varied in the experimental process for each of the materials examined: (1) layer depth, (2) and the nature of the underlying "halfspace". The layer depth was varied between 2 cm and about 40 cm. For the underlying surface, either an aluminum sheet (E = 1 - j3.97 x 107 at 16 GHz) or a flat absorber was used. The flat absorber was composed of 60 x 60-centimeter slabs of 5-centimeter thick 11

flat absorber. These slabs were then arranged to cover the bottom of the sandbox. The dielectric constant of the absorber (E = 1.64 - jO. 15) was measured by comparing the nadir-viewing reflectivity relative to the conducting sheet. The imaginary part was estimated (a lower limit) by considering the minimum amount that the reflection from a conducting sheet under the absorber had been attenuated. Both the radar operation and the movement of the turntable were under computer control. To reduce the variance in the estimate of the mean backscatter, independant measurements were realized by rotating the sandbox. It was found that the signal was sufficiently decorrelated for a rotation of five degrees. Thus, for each target scenario examined (incidence angle, depth, underlying dielectric) 72 independant spatial samples were collected, one every five degrees for one complete revolution of the sandbox. Additional averaging was obtained using frequency averaging, normally using the responses from 21 frequency points equally spaced throughout the 2 GHz bandwidth. The gating capability of the network analyzer was used to isolate the backscatter response of the layer. The footprint of the antenna at 3 meters, even at 60~, the highest incidence angle examined, was sufficiently small (<0.56 m) so that scattering from the wooden sides of the sandbox constituted a negligible interference. Since the primary goal of this investigation is to study volume scattering, great care was taken to make the surface of the layers as smooth as possible to minimize surface backscatter. The smoothing technique employed a sled whose skids rode along two opposite sides of the sandbox, and from which was suspended, spanning the entire width of the sandbox, a metal blade. The metal blade was lowered to the appropriate depth and leveled. The sled was then pulled across the box multiple times. After each pass, excess material was removed manually or additional material added where the level was perceived to be low. In this fashion, an optimally smooth and reproducible surface, with roughness dictated mainly by the particle size of the test material, was achieved. The roughness of the surface was accurately characterized by measurement with a laser profiler. For calibration of the radar, the Isolated Antenna Calibration Technique (IACT) was used [16], which requires the measurement of a sphere (a 6" one in this case) and an arbitrary depolarizing target (we used a length of wire oriented at roughly 45~). This technique, as the name implies, requires a radar system with excellent polarization isolation. The Ku-band system has polarization isolation on the order of 30 dB. 12

5 Application of Hybrid Modeling Approach to Dense Media In this Section, we test the Hybrid Model concept which was proposed in the previous chapter. Specifically we will test whether from measured backscatter data we can retrieve, through an inversion process, parameters which, when used in some form of first-order radiative transfer, provide a comprehensive and polarimetric description of scattering. After evaluating the degree to which such a concept is valid, we will examine the degree to which our measurements of volume scattering agree with theoretical predictions. In general, the inversion process consists of two steps: 1. Assume some form of the unknown matrices 1%, TPbs, and Pbi. A logical starting point is for example, an isotropic model using the three parameters, ie, PI, and P2 from the Rayleigh derivation above. 2. Use an optimization algorithm to minimize the least-squares function given by: No Nd Ne 4 4 [9'jik~lm [5 [.i,k]Im ( 2 i=1 j=1 k=1 1=1 m=l 1,J m where SMU and f7 are the measured and predicted Mueller matrices, respectively, corresponding to all of the different experimental variables-incidence angle, depth, and dielectric of underlying surfaceindexed up to No, Nd, and N,. The formulation of 7 is non-linear, even for the simple case of a scalar extinction, due to the presence of the exponential functions of the extinction. Initially a conjugate gradient procedure was attempted but was found to be unsatisfactory because of the absence of a means for enforcing upper and lower limits on the parameters being estimated. It is, for example, a physical requirement that the parameters corresponding to < SVVS*V >, < ShhS*hh >, and < SvhSv^ > be strictly postitive. We eventually settled upon a very robust algorithm written by the Institute for Systems Research at the University of Maryland [17] called FSQPD (Fortran, Sequential Quadratic Programming, Double precision). The algorithm allows upper and lower bounds to be set for the parameters and also provides for the imposition of linear and non-linear constraints, though none were used in this case. 13

6 Materials Examined The most interesting application of the Hybrid Modeling concept would be to snow. However, in this development stage, the handling requirements of snow for the purpose of a (possible lengthy) set of controlled experiments are too prohibitive. In this initial development stage therefore, we opted to use a "stable" material, for which material characteristics are static and unchanging throughout the layer. The dense medium was a sized silica gravel, henceforth referred to as "8-12 Gravel", based on its being produced in a sieving process as the material retained on Standard Sieves 8-12. Some characteristics of this material are given in Figure which shows the particle size distribution obtained from sieve analyses; also given in the figure are some pertinent physical and electrical properties of the material. The material has a high volume fraction which indicates a maximum packing density of particles. The particle-size distribution of the 8-12 gravel has been rendered as a probability density function to facillitate, as will be discussed below, its comparison with theoretical scattering predictions. The effective dielectric constant of the 8-12 material was measured in two different ways: (i) a coaxial waveguide complex reflection coefficient technique at 9.5 GHz, and (ii) using the Snow Probe [18] at about 1.5 GHz. While the 9.5 GHz result would seem to be most pertinent to the present effort which is conducted at 15-17 GHz, there are questions about the accuracy of the estimate of the imaginary part of the dielectric constant, which is frequency dependent, as measured by the waveguide technique. The issue of determining ~ for the materials is addressed further in Section 9. For now, we show in Figure only the real part of the dielectric constants of the materials, El, which is essentially the same at both of the frequencies at which it was measured. A detailed description of the 8-12 Gravel, provided by the supplier (AGSCO Corp., Wheeling, Il.) is given in Table 6. As shown in Figure b, this material has a relatively narrow particle size distribution, which corresponds to a Gaussian pdf having mean value j = 2.062 mm and variance C2 = 0.056 mm. The particles themselves have a high degree of sphericity. Such factors make it an excellent candidate for testing against theoretical RT models. At Ku-band (X _ 1.875 cm), the associated size parameter, ka, where k is the wavenumber in the particle and a is its radius, is just beyond the Rayleigh limit, for which the criterion has been given as ka < 0.5 [19] (present case: ka =0.61). Although the DMRT model [6] is intended specifically for Rayleigh particles, applications have been demonstrated by its authors for cases having ka as high 14

Table 1: Material characteristics of 8-12 Gravel Typical Physical Properties Typical Chemical Analysis Fusion Point 3135 ~F Silicon dioxide 98.20% Hardness 7 Mhos scale Aluminum oxide 0.72% Grain Shape Spherical Calcium oxide 0.56% Specific gravity 2.65 g/cm3 Iron oxide 0.07% Bulk density 1.57 g/cm3 Magnesium oxide 0.13% pH 6.8-7.2 Manganese oxide 0.03% Sphericity (Krumbein) 0.6-0.7 Potassium oxide 0.10% Roundness (Krumbein) 0.6-0.7 Sodium oxide 0.03% Titanium oxide 0.07% Loss on ignition 0.09% as 1.0 [20]; thus, a meaningful comparison should be permissible with the 8-12 material. One obstacle to theoretical comparisons is the problem of precisely ascertaining the precise complex dielectric constant of the particles themselves at Ku-band. The issue of a theoretical comparision with measurements is taken up in Section 9. 7 Surface Contributions In Section 4 a means by which the top surface of the layer was made as smooth as possible was presented. In order to carefully analyze the scattering behavior of the volume, the effect of the surface must be either (i) very well understood, or (ii) negligible. Since the physics of rough surface scattering is in itself a very challenging problem, and doubly so when integrated with volume scatter, we have attempted to achieve option (ii). The importance of the influence of the surface backscattter on the measurements may be gauged to some degree by the variability of the signal with different target parameters. If, for example, at a given incidence angle, backscatter changes very little with changing depth, or with changes in the underlying "halfspace" material, then it may be suspected that the dominant source of the backscatter is the surface term. As will be shown in the following sections, this was not found to be the case in our measurements. However, even for the case in which the signal is seen to be quite responsive to changes in the target configuration, it is possible that the small surface term may represent a substantial contribution for the cases involving the very lowest backscatter levels, and so introduce errors into an analysis based on volume scattering exclusively. 15

For this reason, an attempt was made to calculate the expected backscatter based on measured parameters of the surface. The surface parameters consist of the effective dielectric constant, given in Figure and the roughness parameters, rms height and correlation length which can be obtained through the measurement of the surface height profile. Accordingly, the surface height profiles were measured with a laser profiler. The RMS height was found to be 0.66 mm, corresponding to a value ks = 0.022. The auto-correlation function was seen to be essentially exponential, with a correlation length of about 2 mm, for which kl = 1.18. A suitable theoretical solution for these parameters is given by the small perturbation method (SPM) (see, for example, [19]). The estimated backscatter as a function of angle is presented in Figure 4. As anticipated, the backscattering level predicted is quite low for this smooth surface. The relative contribution of this surface scattering component may be evaluated with respect to the measurements which are presented next. 8 Interpretation of Results In the following sections we present the results of our controlled experiments and analyze them in the context of the hybrid modeling approach outlined in Section 2. It is our intent to investigate whether the Rayleigh approach described in Section 3 constitutes an appropriate polarimetric model for the test material being considered. Before proceeding further, however, it will first be instructive to state explicitly what consitutes a comprehensive comparison between a polarimetric model and measurements. For an azimuthally symmetric medium, there are normally considered to be five independent quantities contained in the measured Mueller matrix in the backscatter direction. These are the co-pol responses, ISw2 and lShh 2, the cross-pol response, ISvh12 ( = Shv 12 for backscatter) and two parameters which together specify the statistics of the co-polarized phase difference 4hh - v. The pdf governing this random variable has been shown to be [21]: 1 - C2 1 acos(,-;) f(~) 2i[1 -a2cos2(I-.)] /1 2cos2(4)x7 + tan-1 acos(-) ] } (28) which is specified by the two parameters a-the degree of correlation-and p-the co-polarized phase dif 16

ference. The quantities are defined in terms of the elements of the Mueller matrix as, cx- X13+ X14 a=,3 1= tan- X4 (29) V Ml^33 X13 where, 2l = M ' I33 = M22 (30) 2 2 X13 M33+94 14 M34- 43 (31) 44,(1 where 94q are elements of the Mueller matrix. Therefore, a comprehensive comparison will examine the agreement between the model and the data with respect to these five elements. Extensive data was collected for the 8-12 material introduced above, using the experimental procedure described in Section 4. The 8-12 Gravel was examined for ten total layer depths - four over a conductor and six over an absorbing layer (except at 30~, for which it was only two depths over an absorbing layer and five over a conductor) - for incidence angles ranging from 20~ to 60~. There was some indication in the results, however, not predicted by the surface scattering analysis done in the previous section, that at 20~ there was some surface contribution to the backscatter which was not negligible. Therefore the 8-12 Gravel is analyzed using only the data corresponding to 30~, 40~, and 60~. 8.1 Comparison with Rayleigh Model The Rayleigh model derived in Section 3 utilized a scalar extinction, and phase matrix components Pbs and Pbi given by equations (25) and (26). As mentioned, these phase matrix components were entirely specified by just two parameters, which are themselves functions of the elements of the polarizability tensor of the arbitrary Rayleigh particles. This theoretical framework was used along with the results from the measurements in the inversion algorithm described in Section 5. Some selected results of this analysis are shown in Figures 5 and 6. These results are representative of the degree of success which was achieved generally in comparing the Rayleigh model to the measured data. Figure 5 shows the angular variation of the co- and cross-pol responses, examined over both a conductor and an absorber at the depths specified in the figures. In general, the angular variation built into the Rayleigh model does not agree with the observed behavior very 17

well. Figure 6 presents the results with respect to the phase statistics. Only the degree of correlation (ca) is shown at specific layer depths over both absorber and conductor. It is clear that the Rayleigh model is inadequate for explaining the observed behavior for these materials. The lack of success of this model is not too surprising. As mentioned, the individual particles the medium are for the most part pushing the limit of Rayleigh particles for the wavelength being used. When consideration is made of some "effective" aggregate particle, which might comprise many individual particles, the hypothesis becomes even more doubtful. Still it is possible that such a model might be valid for smaller particles and where the volume fraction is not quite so high. The angular variation in the Rayleigh model is controlled by two factors: (1) the Fresnel transmissivity of the surface, and (2) the explicit dependence shown in equation (26). It is possible that a general isotropic model could work if the explicit constraints of the Rayleigh model are relaxed. 8.1.1 General Isotropic Model For the general isotropic model, the quantities of interest take the following forms: * The diagonal elements of the exinction matrix are required to be identical. This is easily understood from the very definition of isotropic media. The propagation characteristics of the medium are insensitive to the polarization of the intensity. Thus the extinction for vertical polarization must be the same as for horizontal. For the non-diagonal elements of the extinction matrix, which for general non-spherical particles is given in [10], we assume the cross-coupling between v and h polarizations is zero for forward scattering, making these off-diagonal elements zero. Thus, for the general isotropic model we are considering, the extinction ice becomes a scalar, ice, as in the Rayleigh model. * The form of!Pbs also resembles that of the Rayleigh model. The isotropic constraint requires that < SSV > =< ShhShh >. Reciprocity requires that Svh = Svh (in the Backscattering Alignment (BSA); in the Forward Scattering Alignment, Svh = -Shv). From arguments of azimuthal symmetry and experimental evidence, the co-pol terms and the cross-pol terms have been found to be statistically uncorrelated; thus all terms of the form < SppSpq > are zero. The resulting form of the matrix, which is 18

independent of incidence angle, is: P2 P3 0 0 P3 P2 0 0 Tbs = (32) 0 0 P4-P3 Ps 0 0 -P5 P4+IP3 * For the Pbi matrix, there is less that can be specified in advance. Reciprocity does not require the Sh and Shv terms to be identical, and the like- and cross-pol terms are not necessarily uncorrelated. For the latter point, we will borrow from the findings of the Rayleigh analysis, which showed these covariance terms to be zero, but relieve generally the other constraints which that analysis placed upon the elements of this matrix. In particular, we will allow the matrix as a whole to be a function of the incidence angle as specified by po (go = cos 0). The structure is, P6 P7 0 0 P8 P9 0 0 Pbi = (33) 0 0 Plo P1i 0 0 P12 P13 The resulting model is identical to that given in Eq. (4), using Pbs and Tbi as given above. There are thirteen parameters to be specified for a complete polarimetric description for a single angle. Of these, the extinction Ke and the matrix Pbs will be common to all angles. A separate bistatic matrix Tbi must be determined for each angle. Thirteen parameters is a considerable space to explore for an algorithm for non-linear optimization. Given the task of optimizing all thirteen parameters at once the algorithm will generally not produce a very satisfactory result, and will tend to arrive at different solutions depending on the initial guess. Fortunately, it is not necessary to optimize all parameters simultaneously. The following option exists for a much more limited optimization process. First, it is recognized that, if a "quadrant" is defined by a 2 x 2 submatrix, then proceeding clockwise from the top left of the 4 x 4 matrices in this formulation, the first and third quadrants are decoupled from 19

one another, since the second and fourth quadrants are identically zero. The zero status of these latter two quadrants holds for all of the matrices involved in the first-order solution, including the transmissivity and reflectivity matrices. Thus, it is possible to optimize initially only the seven parameters affecting the coand cross-pol responses. After these have been found, the remaining six may be found, using the value for Ke obtained in the initial process. Additional iterations may be performed to improve the overall result. In general, it is found that, due to a model which is only an approximation and data which contains errors, no absolute convergence is observed from an iterative process, and degradation of the perceived "goodness" of the solution occurs with many iterations. A complete comparison of the isotropic modeling approach with the measured data is shown in Figures 7-9 for the 8-12 Gravel, for incidence angles 300, 400, and 600. All of the available data was used in the inversion process suggested by Eq. (27 to find the optimal parameters for the first-order polarimetric model given by Eq. (4). Each of the figures constitutes a comprehensive comparison of the five aforementioned independent elements of the measured Mueller matrix with the model predictions for the angle specified. The agreement with the 8-12 Gravel is seen to be generally very good. The co- and cross-polarized response for layers over a conducting surface ((a) in the figures) is very well modeled. There is some small disagreement for the case of the co-pol response over an absorber layer, at all three angles. The response is somewhat underpredicted (by 1-2 dB) at 300, and over predicted by about this amount at 60~. This is an interesting and important case to examine, since the bistatic component of the phase function, Tbj, which can be individually tailored for each incidence angle in the inversion process, is not involved. Since the backscatter component, fPbs, is a constant quantity, the angular variation is governed solely by the Fresnel transmissivity of the surface. Angular behavior which is in agreement with this Fresnel variation, for a "halfspace", which our layers over absorber approximate, is a necessary requirement for scattering which is truly first-order and isotropic. The behavior of the phase statistics of the measurements appears to be a very simple function of all of the parameters varied (angle, depth, underlying "halfspace"), and is modeled fairly well by this first-order hybrid model. The parameters 4, the co-polarized phase difference is, in particular, explained very well by the is modeled quite well at all three angles. 20

9 Comparison of Results with Theory In Section 6 we stated that the characteristics of the 8-12 Gravel made it especially suitable for comparison with theory. One impediment, however, to such a comparison is the task of determining the precise dielectric constant of the particles themselves, Ep, at the frequency of interest, in this case, our center frequency of 16 GHz. One approach by which this might be accomplished is to measure by some means the effective dielectric constant, eff, of the medium and then invert a dielectric mixing formula like the Polder-Van Santen mixing model [19] to get the dielectric of the inclusions (particles). The fact that the particles have a high degree of sphericity removes the ambiguity associated with assigning the three shape factors for that model. If, however, Eeff is measured at the test frequency (16 GHz was the center frequency), then the imaginary part, ~eff9 may comprise the effects of the scattering losses as well as the dielectric conducting losses depending on the measurement technique. If the scattering losses are included in the measurement of Eeff, inversion of a dielectric mixing model will not lead to the proper result for the particle dielectric Ep. A complex reflection coefficient waveguide technique might be employed at this frequency (16 GHz) in which (it may be argued) scattering losses are prevented due to the boundaries (i.e. modal requirements) of the guide; this point however, is at this time still something of a research question, and in any event, the technique is prone to error when the loss tangent of the test material is very small. A more appropriate procedure is to measure E" at low frequency, where scattering is negligible independent of the measurement technique used. The variation of the dielectric constant of dry rocks with frequency in the microwave region has been studied [22]. The real part p' is essentially constant from 1 to 16 GHz. The imaginary part ep has been found to decrease with frequency, in a manner dependent on the particular rock class. The material in the 8-12 Gravel appears to be either a sedimentary, plutonic, or volcanic silicate. The frequency dependance of ep for each one of these classes is given in [22]. In this case, Eeff for the 8-12 Gravel was measured at about 1.5 GHz using the Snow Probe (described in [18]) and found to be Leff = 2.17 - j0.007. The Polder Van Santen mixing formula for spherical particles (shape factors: AI = A2 = A3 = 1/3) is, em = Eh + 3viEm e -,E E( +}-2EMi 21

where Em, Eh, and ei are the effective dielectric constants of the medium, the host (air), and the inclusions (or particles), respectively, and vi is the volume fraction. This equation may be easily solved for ei (= Ep) yielding 3.11 - j1.38 x 10-2. This value of E is already considerably lower than the average for all of the silicate classes mentioned above (see Figure 13 in [22]). Therefore, instead of applying the frequency dependance associated with the much higher values of the quantity E'p measured for the various silicate classes in [22], we will simply use the low frequency value of Ep as our estimate at the frequency of interest, 16 GHz. This result completes the information required to examine the correspondence between the measurements of this material and the predictions of certain discrete-particle-based theories. A pertinent quantity for comparison is the effective propagation constant in the medium: K = K'- jK" from which may be obtained the real and imaginary parts of the effective index of refraction: K' K" n = n"= (34) ko 'o with the extinction specified as ie = 2K". The theories we compare with are the Effective Field Approximation (EFA), the Quasi-Crystalline Approximation (QCA), and the Quasi-Crystalline Approximation with Coherent Potential (QCA-CP). We present here a very brief description of each technique. In conventional radiative transfer (CRT) the extinction is obtained using EFA (also known as Foldy's approximation), which can be considered a special case of QCA, [14] in that the particles are considered as acting completely independent of one another with zero correlation. The solution for the effective propagation constant is formally given by [14], K = [k2 - 4nno < Spq(Os, s; i, Xi) >]1/2 where Spq is the complex far-field scattering amplitude, k is the wavenumber of the surrounding medium and no is the number density of the scatterers. The loss the wave experiences is due to the total extinction cross-section of a particle multiplied by no. 22

Table 2: Theoretical predictions based on 8-12 Gravel characteristics. Source n' K Meas. 1.47 2.35 EFA 1.35 3.09/3.02* QCA 1.43 0.29 QCA-CP~ 1.47 1.69 *Background = air / Background = Eeff 'Intended for Rayleigh region. A detailed description of QCA theory is outside the scope of this presentation. A detailed derivation may be found in [14]. In general, QCA takes into account interactions between particles, using exact wave transformations, a "T-matrix" approach, for computing this interparticle interaction rather than far-field phase functions. The particle correlation is specified by the pair-distribution function which describes the conditional probability of a particle location relative to another particle's position. The Percus-Yevick approximation is most often used to derive the pair-distribution function for a medium. This approximation assumes non-interpenetrability of particles and zero forces between the particles. This assumption, which has, in some specific experiments [23], been found to agree with experimental observations, removes the requirement of characterizing by more direct means the exact configuration of particles in a dense medium. In more recent studies [24] this assumption has been shown to break down outside of its intended domain of validity, which comprises mainly liquids or gases. QCA-CP constitutes an improvement upon QCA in that it provides for energy conservation in the formulation. Essentially it amounts to employing the effective wavenumber K in the Green's functions of the QCA formulation instead of the background wavenumber k. The formulation is valid for Rayleigh particles only. The results produced by each one of these methods are summarized in Table 2. Shown are estimates for n', the real part of the effective index of refraction and Ke. Shown for comparison, as "Meas.", are the Snow Probe measured n' and Ke obtained from the isotropic Hybrid model analysis of the 8-12 Gravel. The correlated particle treatments represented by QCA and QCA-CP produce estimates of n' which are in close agreement with the measured value. EFA gives a considerably lower estimate. While the QCA result for Ke is very low, the QCA-CP and EFA estimates bracket, below and above, respectively, the experimentally determined value for Ke. This bracketing scenario resembles previous findings [23] for glass spheres in styrofoam 23

at much lower volume fractions (< 10%). Figure 10 shows a comparison between measurements and scattering computations from a numerical CRT model which employs the EFA theory, using Mie calculations for the phase function and extinction and the discrete ordinate solution method [25]. The inputs to the model consisted of the physical characteristics of the 8-12 gravel as they have been described in the preceding sections, including particle size distribution, volume fraction, effective dielectric constant of the medium and dielectric constant of the individual particles. Though the trends with respect to polarization and depth are similar, the CRT model overestimates the scattering level by typically 3-4 dB. A comparison between measurements and predictions from a Dense Medium radiative transfer (DMRT, described in Section 2) model are shown in Figure 11. The theoretical predictions are seen to be extremely low-on the order of 15-20 dB down-relative to the measured results. This very serious disagreement can be explained in terms of the behavior of the scattering albedo, ca, which scales the phase matrix (as described in Section 2) and hence the scattering production of a collection of scatters. The value of co computed (from [8]) for the DMRT model is similar to the CRT computation for very low volume fractions (eg. <5%) but diminishes rapidly with increasing volume fraction. For this scenario, with a volume fraction of 0.63, the scattering albedo computed from DMRT is about 250 times smaller than that computed by CRT. Since the values of extinction given in Table 2 for EFA and QCA-CP are fairly similar, this means that the difference in the respective albedos pertains to the scattering (recall, Co, = Ky/KY). The value for the extinction computed from QCA-CP given in Table 2 is based almost entirely on the absorptive losses in the particles. That the QCA result is smaller than this value is apparently symptomatic of that technique's failure to obey energy conservation. For a final comparison, we note that in the isotropic Hybrid model, finding Tib independently for various incidence angles amounts to mapping out some portion of the scattering phase function of the particle (refer, for example, to Figure 1). Figure 12 shows these experimentally derived phase functions. Shown for comparison are Mie calculations made for the cases of (a) a particle having the same characteristics as the 8-12 Gravel particles, and (b) a particle having the same dielectric but three times as large. All results are normalized to the backscatter level. The obvious implication, consistent with the assumptions of the Hybrid model concept, is that multiple particles combine to form an effective scattering element which is amenable to a first-order RT approach. 24

10 Summary and Conclusions This paper has described a hybrid experimental / theoretical technique by which dense medium scattering may be modeled. Specifically, we have demonstrated that the complete polarimetric backscatter response from a very dense medium consisting of discrete particles can be adequately explained in terms of a vector first-order radiative transfer model. Very simple forms of the phase matrix components involved in the first-order formulation have been derived using Rayleigh theory. While these specific formulations were not observed to correctly model the experimentally observed results in the present case, it is possible they may still be valid if applied to a somewhat less dense medium with smaller particles. The dense medium backscatter measurements were compared against predictions from conventional and dense media radiative transfer. While the CRT results were reasonably close, the scattering predicted by DMRT for this high volume fraction material (vi = 0.63) was 15-20 dB below the measured levels. Finally, application of the hybrid technique to direct measurements though a numerical optimization algorithm also yields a partial mapping of the scattering phase function of the effective volume scattering element in the medium. An analysis has been presented which suggests that this element is considerably larger than the physical particles of which the medium is comprised. 25

References [1] A. Nashashibi and K. Sarabandi, "Propagation in dense random media: A novel measurement technique and experimental results," IEEE Trans. Antennas Propagat., 1995. submitted for publication. [2] J. R. Kendra, K. Sarabandi, and F. T. Ulaby, "Radar measurements of snow: Experiment and analysis," IEEE Trans. Geosci. Rem. Sens., 1995. submitted for publication. [3] F. Vallese and J. A. Kong, "Correlation function studies for snow and ice," J. Appl. Phys., vol. 52, Aug. 1981. [4] G. Koh and R. E. Davis, "Effect of snow stereology on millimeter wave extinction," in Proc. IGARSS '94, p. 1929, 1994. Vol. IV. [5] P. R. Siqueira, K. Sarabandi, and F T. Ulaby, "Numerical simulation of scatterer positions in a very dense medium with an application to the two-dimensional born approximation," Radio Science, 1995. accepted for publication. [6] L. Tsang and A. Ishimaru, "Radiative wave equations for vector electromagnetic propagation in dense nontenuous media," J. Electromagnetic Waves and Applications, vol. 1, no. 1, pp. 59-72, 1987. [7] F T. Ulaby and C. Elachi, eds., Radar Polarimetryfor Geoscience Applications. Dedham, MA: Artech House, Inc., 1990. [8] B. Wen, L. Tsang, D. P. Winebrenner, and A. Ishimaru, "Dense medium radiative transfer theory: comparison with experiment and application to microwave remote sensing and polarimetry," IEEE Trans. Geosci. Rem. Sens., vol. 28, Jan. 1990. [9] L. Tsang and J. A. Kong, "Scattering of electromagnetic waves from a dense medium consisting of correlated mie scatterers with size distributions and applications to dry snow," Journal of Electromagnetic Waves and Applications, vol. 6, no. 3, pp. 265-296, 1992. [10] A. Ishimaru and R. L.-T. Cheung, "Multiple scattering effects on wave propagation due to rain," Ann. Telecommunication, vol. 35, pp. 373-379, 1980. 26

[11] M. W. Whitt, Microwave scattering from periodic row-structured vegetation. PhD thesis, The University of Michigan, Ann Arbor, MI, 1991. [12] R. E. Kleinman and T. B. A. Senior, "Rayleigh scattering," in Low and High Frequency Asymtotics (V. K. Vardan and V. V. Varadan, eds.), ch. 1, Amsterdam: North-Holland, 1986. [13] J. R. Kendra, Microwave Remote Sensing of Snow: An Empirical / Theoretical Scattering Model for Dense Random Media. PhD thesis, The University of Michigan, 1995. [14] L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing. New York, NY: Wiley Interscience, 1985. [15] M. A. Tassoudji, K. Sarabandi, and F T. Ulaby, "Design consideration and implementation of the lcx polarimetric scatterometer (POLARSCAT)," Tech. Rep. Rep. 022486-T-2, The University of Michigan, Radiation Lab, 1989. [16] K. Sarabandi, F T. Ulaby, and M. A. Tassoudji, "Calibration of polarimetric radar systems with good polarization isolation," IEEE Trans. Geosci. Rem. Sens., vol. 28, no. 1, 1990. [17] J. L. Zhou and A. L. Tits, "User's guide for fsqp version 3.3b: A fortran code for solving constrained nonlinear (minimax) optimization problems, generating iterates satisfying all inequality and linear constraints," Tech. Rep. TR-92-107r3, Systems Research Center, The University of Maryland, September 1993. [18] J. R. Kendra, F. T. Ulaby, and K. Sarabandi, "Snow probe for in situ determination of wetness and density," IEEE Trans. Geosci. Rem. Sens., vol. 32, no. 6, 1994. [19] F T. Ulaby, R. K. Moore, and A. K. Fung, Microwave Remote Sensing, Active and Passive, vol. 1-3. Dedham, MA: Artech House, Inc., 1986. [20] R. West, L. Tsang, and D. P. Winebrenner, "Dense medium radiative transfer theory for two scattering layers with a rayleigh distribution of particle sizes," IEEE Trans. Geosci. Rem. Sens., vol. 31, pp. 426 -437, Mar 1993. [21] K. Sarabandi, "Derivation of phase statistics from the mueller matrix," Radio Science, vol. 27, no. 5, pp. 553-560, 1992. 27

[22] F. T. Ulaby, T. H. Bengal, M. C. Dobson, J. R. East, J. B. Garvin, and D. L. Evans, "Microwave dielectric properties of dry rocks," IEEE Trans. Geosci. Rem. Sens., vol. 28, no. 3, 1990. [23] C. E. Mandt, "Microwave propagation and scattering in a dense distribution of spherical particles," Master's thesis, Univ. of Washington, Seattle, 1987. [24] P. R. Siqueira and K. Sarabandi, "Numerical evaluation of the two-dimensional quasi-crystalline approximation," IEEE Trans. Antennas Propagat., 1995. accepted for publication. [25] Y. Ulaby, T. F Haddock, and R. D. DeRoo, "Millimeter-wave radar scattering from snow: 1. Radiative transfer model," Radio Science, vol. 26, no. 2, pp. 329-341, 1991. 28

List of Figures 1 First-order volume scattering mechanisms in a layer of scatterers................ 30 2 Photograph of experimental setup. Radar mounted on a mechanical lift overlooking a sandbox filled with the test material. (Additional equipment is part of an unrelated bistatic facility.) 30 3 Characteristics of (a) Common sand, and (b) 8-12 Gravel: vi denotes volume fraction. Note that distribution in (b) has been rendered as a probability density function........... 31 4 Predictions of surface scattering for 8-12 Gravel, generated using the Small Perturbation M ethod............................................. 31 5 Application of Rayleigh model to measured data. Comparison for co- and cross-polarized response, over both absorber and conductor for 8-12 gravel.................. 32 6 Application of Rayleigh model to measured data. Comparison for degree of correlation (a). Layers over absorber and conductor are 14.3 and 16.5 cm respectively................. 33 7 Results of application of Isotropic RT model to measured data: 8-12 Gravel, 30~ incidence.. 33 8 Results of application of Isotropic RT model to measured data: 8-12 Gravel, 40~ incidence.. 34 9 Results of application of Isotropic RT model to measured data: 8-12 Gravel, 600 incidence.. 35 10 Comparison of measurements of 8-12 Gravel with CRT simulations. Discrete marks are CRT simulations; continuous lines are from non-isotropic, first-order RT analysis of measured data which very closely approximates the actual data...................... 36 11 Comparison of measurements of 8-12 Gravel with DMRT simulations. Discrete marks are DMRT simulations; continuous lines are from non-isotropic, first-order RT analysis of measured data which very closely approximates the actual data................... 36 12 A portion of the scattering phase function obtained via isotropic RT analysis of 8-12 Gravel compared with Mie phase function for a particle (a) identical to an 8-12 Gravel particle, (b) three times larger than an 8-12 Gravel particle.......................... 37 29

Figure 1: First-order volume scattering mechanisms in a layer of scatterers. Figure 2: Photograph of experimental setup. Radar mounted on a mechanical lift overlooking a sandbox filled with the test material. (Additional equipment is part of an unrelated bistatic facility.) 30

2r >1 -l c.0 co a. =- 2.06 mm.5 - (c =0.24 mm 1.5 -1 0 1 I 8-12 Gravel ~'= 2.17 v.= 0.63 11 2 3 Diameter (mm) I I I 4 5 6 J 7 Figure 3: Characteristics of (a) Common sand, and (b) 8-12 Gravel: vi denotes volume fraction. Note that distribution in (b) has been rendered as a probability density function. 0.0 -10.0 -20.0 -30.0 -40.0 ~~0 0 -50.0 L' 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Incidence Angle (deg.) Figure 4: Predictions of surface scattering for 8-12 Gravel, generated using the Small Perturbation Method 31

0.0 -5.0 -10.0 -15.0 m ~0 -20.0 -25.0 - -30.0 K 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (deg.) (a) 14.3 cm over absorber. 0.0 -5.0 -10.0 -15.0 "10 03 -20.0 - -25.0 - -30.0 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (deg.) (b) 16.5 cm over conductor. Figure 5: Application of Rayleigh model to measured data. Comparison for co- and cross-polarized response, over both absorber and conductor for 8-12 gravel. 32

1.00................ 0.80 - - CcA o 0.40 — 0- OVERABS J..0. - OVER PEC 0.20 - - 0.00 I...I 20.0 30.0 40.0 50.0 60.0 70.0 Incidence Angle (deg.) Figure 6: Application of Rayleigh model to measured data. Comparison for degree of correlation (a). Layers over absorber and conductor are 14.3 and 16.5 cm respectively. -— o --- VV -10.0 -:-::X —:-3-:1g: -.-0.0 - | | _ -5.0 -5.0 - -- -- tHH 0 -20.0 0 -1 0.0 - 1,...,....... 20,,,. -... - -30.0- -30-.0 -..-35.0... —.......... - 0.0 10.0 20.0 30.0 40.0 50.0 0.0 10.0 20.0 30.0 40.0 50.0 Layer Depth (cm) Layer Depth (cm) (a) Co- and Cross over conductor. (b) Co- and Cross over absorber. 1.00 i I i 200 0.90 - 195 1 —0 — OVER ABS 0.80 - - 190 -000,,,, O....;.,,C, 0.70 - -................ 185 - -. -- -------- --— ~~ - 0.60 180 - --------------------------------- 0.50 175 —. --- —-—.-. ----- ------------- o 0.40 -30 - 0.30 - 165 - 2 — 0 --- Overabs. 1 —.. --- Overpec 0.10 - 155 0.00....150....... I..... 0.0 10.0 20.0 30.0 40.0 50.0 60.0 0.0 10.0 20.0 30.0 40.0 50.0 60.0 Layer Depth (cm) Layer Depth (cm) (c) Degree of Correlation. (d) Co-Polarized Phase Difference. Figure 7: Results of application of Isotropic RT model to measured data: 8-12 Gravel, 30~ incidence. 33

0.0. *... i... *...... *.... * *. ' * w... * | * *. * |,..., * * *. 0.0 I I I -5.0 - -10.0 -10.0-, - * - -' ' '-*- -.-:-. —:.:.-:-3- '-: —::'.:::::'-_:.:-:: |-10.0 a '. " 25 -,'."'. ->~A A -25.0 /.o- vv -30.0:0.90- -- -- VV / — 3 — HH AB /0.0 - - -- HH -30.0 - --— B, --- VH H —.. ---- VH -35.0 -40.0.. 0.0 10.0 20.0 30.0 40.0 50.0 0.0 10.0 20.0 30.0 40.0 50.0 Layer Depth (cm) Layer Depth (cm) (a) Co- and Cross over conductor. (b) Co- and Cross over absorber. 1.F00i i. I d i.. 200... 0.90 - 195 -— o --- OVER ABS 0.80 - 190 - o o '..-..B. OVERPEC 0.70 -..'''" 18-8 ---. —............. 185 - ----------------------------— 0 0.60 - 180... --- —----—, ----—...' — ---- o 0.50 - 175........o 0 -{7.............................. o 0.40 - 170 - 0.30 - 165 - ' — 0 --- Overabs. *' 0.20 a- 160 -..... Over pec 6 0.10 155 0.00 150 I.... 0.0 10.0 20.0 30.0 40.0 50.0 60.0 0.0 10.0 20.0 30.0 40.0 50.0 60.0 Layer Depth (cm) Layer Depth (cm) (c) Degree of Correlation. (d) Co-Polarized Phase Difference. Figure 8: Results of application of Isotropic RT model to measured data: 8-12 Gravel, 40 incidence. 34

-10.0 -20.0 0o ~G -30.0 -40.0 - 0.0 50.0 0 0 0 0 0 ee 00 0 ao Layer Depth (cm) (a) Co- and Cross over conductor. 1.00......... 0.90 0.80 0.70 0 o 0.60 0 ro o ~~~....0.50............... a... 0.40 0.30 ------ Over abs. 0.20 ~-.... Over pec 0.10 0.00.... 0.0 10.0 20.0 30.0 40.0 50.0 60.0 Layer Depth (cm) (c) Degree of Correlation. L.P I4) lu 4) Layer Depth (cm) (b) Co- and Cross over absorber. 200............. 195 -— o --- OVERABS 190 o ~. --- OVER PEC 185 - 180 -------- -- ----------------. 175.170 0 -165 160 155 150.. 0.0 10.0 20.0 30.0 40.0 50.0 60.0 Layer Depth (cm) (d) Co-Polarized Phase Difference. Figure 9: Results of application of Isotropic RT model to measured data: 8-12 Gravel, 60~ incidence. 35

0.0 o -10.0.......gI ~"i....i.... -- ---- VV 88g @ - *- - HH o'............. /I -10.0 0 eP ~o -20.0 -30.0 0. 0....................... 10.0 20.0 30.0 40.0 50.0 Layer Depth (cm) (a) 30~ over conductor. 0.0 10.0 20.0 30.0 40.0 50.0 Layer Depth (cm) (b) 300 over absorber. 0.0 -10.0 -20.0 — 9 --- VV - --- HH -.. ---- VH -...................... /..,. ^,.- ^- - - -.. -.. -10.0 I30 -20.0 P.............. A A? A] l - -..... --- —------ --- - -.. —. - /' — <a — HH / —..- VH / -30.0 -30.0 L 0. 0 -... I.... I -... I.... ]li% I%................ 10.0 20.0 30.0 40.0 50.0 0.0 10.0 20.0 30.0 40.0 50.0 Layer Depth (cm) (c) 60~ over conductor. Layer Depth (cm) (d) 60~ over absorber. Figure 10: Comparison of measurements of 8-12 Gravel with CRT simulations. Discrete marks are CRT simulations; continuous lines are from non-isotropic, first-order RT analysis of measured data which very closely approximates the actual data. -10.0 -20.0 ~., ------— " /..,...... vv -o!- HH I n / ti -30.0 P -40.0 L0.0 -40.0 L O.0... I.... I.... I.... I.. 10.0 20.0 30.0 40.0 50.0 Layer Depth (cm) (a) 30~ over conductor. 10.0 20.0 30.0 40.0 50.0 Layer Depth (cm) (b) 30~ over absorber. Figure 1 1: Comparison of measurements of 8-12 Gravel with DMRT simulations. Discrete marks are DMRT simulations; continuous lines are from non-isotropic, first-order RT analysis of measured data which very closely approximates the actual data. 36

ce - O).' 0 C 3 o P1 c. 4.00 3.00 I I I I "l I I -0 — HH (MEAS).......... HH (MIE)..A-.- MA,. A -I A —A --- \ VV (MEA) - VV (MIE) 2.00 \'4 i......................,....o ~.b~............ "' "...,., --... \\ -^^^^^^^^ 1.00 0.00 L * i 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Bistatic Scattering Angle (deg.) (a) 4.00., —4 U).r~ 3.-1 rc I.. GO 3.00 2.00 1.00 0.00 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 Bistatic Scattering Angle (deg.) (a) Figure 12: A portion of the scattering phase function obtained via isotropic RT analysis of 8-12 Gravel compared with Mie phase function for a particle (a) identical to an 8-12 Gravel particle, (b) three times larger than an 8-12 Gravel particle. 37