# Application of Passive and Active Microwave Remote Sensing for Snow Water Equivalent Estimation

код для вставкиСкачатьApplication of Passive and Active Microwave Remote Sensing for Snow Water Equivalent Estimation Dissertation Presented in Partial Fulfillment of the Requirements for the Degree Doctor of Philosophy in the Graduate School of The Ohio State University By Jinmei Pan, M.S. Graduate Program in Geodetic Science The Ohio State University 2017 Dissertation Committee: Michael T. Durand, Advisor Che-Kwan Shum Ian M. Howat Joel T. Johnson Barbara E. Wyslouzi ProQuest Number: 10702692 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. ProQuest 10702692 Published by ProQuest LLC (2017 ). Copyright of the Dissertation is held by the Author. All rights reserved. This work is protected against unauthorized copying under Title 17, United States Code Microform Edition © ProQuest LLC. ProQuest LLC. 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, MI 48106 - 1346 Copyrighted by Jinmei Pan Abstract Snow accumulation on the ground changes the energy balance between the land and the atmosphere, and stores winter precipitation for water supplies in many parts of the world. In practice, the snow water equivalent (SWE), defined as the equivalent depth of liquid water when snow completely melts, is difficult to map in cold regions except via remote sensing techniques. The microwave remote sensing (MWRS) has been used for SWE estimation since the 1980s based on the interactions of microwave radiation with snow crystals. In this study, physically based radiative transfer (RT) models and the Bayesianbased Markov Chain Monte Carlo (MCMC) approach were applied to develop a highaccuracy SWE retrieval algorithm. The models and the algorithms were tested using ground-based snowpit and microwave measurements. Two widely-used snow RT models were fully-compared in the aspects of snow micro-structure assumptions, volume scattering theories and the RT equation resolution. The Microwave Emission Model of Layered Snowpacks (MEMLS) based on the Improved Born Approximation (IBA) was shown to be an adequate observation model to estimate SWE using the multi-frequency brightness temperature (TB) at 10.65 to 90 GHz. The prior information is from a set of globally-available datasets, and the performance is tested for local prior information derived from historical ground measurements. The retrieval algorithm was later adapted for active microwave SWE retrieval using the backscattering coefficient at 10.2 to 16.7 GHz. Results showed that MEMLS-IBA can simulate the measured microwave signals ii with a 10-K accuracy for TB and a 1-dB accuracy for the backscattering coefficient. The passive microwave retrieval algorithm achieved an accuracy of 30-mm for shallow snow, with two-layer snow properties estimated. However, the active microwave retrieval algorithm reproduced similar accuracy only in the synthetic experiment using 1-layer snow property estimates. Future improvement requires a better active microwave observation model. iii To my family, my advisors and my friends. iv Acknowledgments I want to sincerely acknowledge all the individuals who dedicated to the Nordic Snow Radar (NoSREx) Experiements, the Cold Land Processes Field Experiments (CLPX), and the Canadian CoReH2O Snow and Ice Experiments, and collected enormous valuable datasets to be used by the community. I want to express my full gratitude to my advisor and labmates for the assistance they kindly provided to support my Ph.D. research. Thanks to my family and my friends, for the happiness they brought to my life. This work is supported by the funds from the China Scholarship Council (2012-2016), the OSU Presidential Fellowship (2016-2017) and the NASA New Investigator Program (NNX13AB63G). v Vita 2005................................................................Chongqing Foreign Languages School 2009................................................................B.S. Geographical Information System, Beijing Normal University 2012................................................................M.S. Cartography and Geographical Information System, Beijing Normal University 2012 to present ..............................................School of Earth Science, The Ohio State University Publications Jinmei Pan, Michael Duarnd, Benjamin Vander-Jagt & Desheng Liu. (2017). Application of a Markov Chain Monte Carlo Algorithm for Snow Water Equivalent Retrieval From Passive Microwave Measurements. Remote Sensing of Environment, 192: 150-165. Jinmei Pan, Michael Durand, Melody Sandells, Juha Leymmetyinen, Edward Kim, Jouni Pulliainen, Anna Kontu & Chris Derksen (2016). Differences Between the HUT Snow vi Emission Model and MEMLS and Their Effects on Brightness Temperature Simulation. IEEE Transactions on Geoscience and Remote Sensing, 54(4): 2001-2019. Jiuquan Cheng, Lingmei Jiang, Lixin Zhang & Jinmei Pan (2015). Microwave Emission Characteristics of Snow Melting/Refreezing Cycles in North China. Remote Sensing Information (in Chinese), 3: 65-73. Xiyu Chen, Lingmei Jiang, Juntao Yang & Jinmei Pan (2014). Validation of Ice Mapping System Snow Cover over Southern China Based on Land Enhanced Thematic Mapper Plus Imagery. Journal of Applied Remote Sensing (in Chinese), 8(1): 084680. Juntao Yang, Lingmei Jiang, Jinmei Pan & Lixin Zhang (2013). Study of Monitoring Snow Cover Using GOES-Geostationary Satellite and AMSR-E. Remote Sensing Technology and Application (in Chinese), 28(5): 920-927. Jinmei Pan, Lingmei Jiang & Lixin Zhang (2012). Comparison of the Multi-layer HUT Snow Emission Model with Observations of Wet Snowpacks. Hydrological Processes, 28(3): 1071–1083. Jinmei Pan, Lixin Zhang, Haoran Wu, Junping Zhuang, Shaojie Zhao, Tianjie Zhao & Lingmei Jiang (2012). Effect of Soil Organic Substance on Soil Dielectric Constant. Journal of Remote Sensing (in Chinese), 16(1): 1-12. vii Yu Liu, Lingmei Jiang, Jiancheng Shi, Shenglei Zhang, Jinmei Pan & Pei Wang (2011). Validation and Sensitivity Analysis of the Snow Thermal Model (SNTHERM) at Binggou basin, Gansu. Journal of Remote Sensing (in Chinese), 2011, 5(4): 1993-2002. Fields of Study Major Field: Geodetic Science viii Table of Contents Abstract ............................................................................................................................... ii Acknowledgments............................................................................................................... v Vita..................................................................................................................................... vi Publications ........................................................................................................................ vi Fields of Study ................................................................................................................. viii Table of Contents ............................................................................................................... ix List of Tables .................................................................................................................... xii List of Figures .................................................................................................................. xiv Chapter 1: Introduction ....................................................................................................... 1 Section 1.1: Snow Radiative Transfer Models ................................................................ 4 Section 1.2: Passive SWE Retrieval Algorithms ............................................................ 8 Section 1.3: Active SWE Retrieval Algorithm ............................................................. 11 Chapter 2: Snow Radiative Transfer Model Comparison ................................................. 14 Section 2.1: Differences in Model Theory .................................................................... 14 One-Flux Versus Two-Flux Theories ........................................................................ 15 ix Scattering Coefficient ................................................................................................ 21 Trapped Radiation ..................................................................................................... 24 Section 2.2: Experiments to Test Three Model Theoretical Differences ...................... 25 Experiment Design .................................................................................................... 26 Simulated TB Differences from Model Theories ....................................................... 28 Section 2.3: Comparison of Model Simulations with Passive TB Measurements ......... 39 Section 2.4: Comparison of Model Simulations with Active Backscattering Measurements................................................................................................................ 54 MEMLS3&a Introduction ......................................................................................... 54 MEMLS3&a Simulation Results ............................................................................... 57 Discussion.................................................................................................................. 63 Chapter 3: SWE Retrieval Using Passive Microwave Measurements.............................. 71 Section 3.1: Estimated Parameters for Passive SWE Retrieval .................................... 72 Section 3.2: The BASE-PM Retrieval Algorithm Theory and Implementation ........... 74 Section 3.3: Priors for SWE Retrieval .......................................................................... 77 Section 3.4: Passive Microwave Retrieval Results ....................................................... 84 Retrieved SWE Using Real TB Measurements .......................................................... 85 Retrieved SWE Using Synthetic TB .......................................................................... 89 Markov Chains of Estimated Variables ..................................................................... 91 x Section 3.5: Discussions................................................................................................ 98 Analysis of Response of Microwave TB to SWE ...................................................... 98 The Influence of The Number of Snow Layers ....................................................... 101 Chapter 4: SWE Retrieval Using Active Microwave Measurements ............................. 108 Section 4.1: Active Microwave Retrieval Algorithm Setup ....................................... 108 Section 4.2: Retrieval Results Using SnowScat Measurements ................................. 109 Section 4.3: Retrieval Results Using Synthetic Backscattering Coefficient ............... 120 Section 4.4: Discussion ............................................................................................... 126 Chapter 5: Conclusions ................................................................................................... 130 Bibliography ................................................................................................................... 133 Appendix A: Solving A0 and B0 For 1-flux RT Equation .............................................. 143 Appendix B: Solving RT Equations For Multi-Layer Snowpack................................... 144 xi List of Tables Table 1. The scattering coefficients calculated by the Improved Born Approximation at 37 GHz for three sizes of snow particles .......................................................................... 26 Table 2. Difference of simulated TB above the snow surface for 40-cm, 1-meter and 2meter slant snow thickness with medium snow grain size at 37 GHz. ............................. 39 Table 3. Details of TB measurements at all sites. .............................................................. 41 Table 4. Mean biases of the simulated compared to the observed TB above the snow surface at all sites .............................................................................................................. 52 Table 5. Rms errors of simulated compared to the observed TB above the snow surface at all sites .............................................................................................................................. 53 Table 6. Mean bias and rms error of the MEMLS3&a simulations for the Sodankyla snowpits. ........................................................................................................................... 63 Table 7. Mean bias and rms error of MEMLS3&a simulations for snowpits not influenced by part-frozen soil conditions. .......................................................................................... 69 Table 8. Statistics of SWE and the mass-weighted average density, temperature and micro-structure parameters of the measured snowpits...................................................... 79 Table 9. Monthly-averaged SWE from the VIC model predictions and the generic priors of snow parameters used for retrieval ............................................................................... 81 Table 10. Error of MCMC-retrieved SWE using observed TB ......................................... 89 xii Table 11. Error of MCMC-retrieved SWE using synthetic TB ......................................... 91 Table 12. The estimated snow and soil parameters for the example Sodankyla snowpit observed on March 2, 2010 when generic SWE and stratigraphy priors were used. The “true” values of the snow variables were relayered from the original 6-layer measurement. .................................................................................................................... 97 xiii List of Figures Figure 1. A homogenous snowpack of thickness, d, where z=0 is at the snow surface, and z=-d is at the bottom of snow. ........................................................................................... 18 Figure 2. Four TB quantities, Aj, Bj, Cj and Dj, inside the j-th snow layer defined in Wiesmann & Matzler (1999a). The z-axis is along the vertically upward direction. zj-1 and zj indicate the lower and upper boundary of the j-th snow layer. In the figure, sgs is the ground-snow boundary reflectivity; sas is the air-snow boundary reflectivity; Tground is the upward ground TB; and, Tsky is the downward sky TB. ..................................................... 21 Figure 3. Upward TB profiles of the 1-flux and 2-flux models using the same absorption coefficient and scattering coefficient. The value of scattering coefficient refers to the Improved Born Approximation (IBA) for snow with medium grain size (pec=0.18 mm) at 37 GHz, and the absorption coefficient refers to snow with a density of 200 kg m-3 and a physical temperature of -5 °C at 37 GHz. The red thick line is the 1-flux model simulation with q=0.96. The red thick dash line is the 1-flux model simulation with q=0.5. The blue square line is the 2-flux model simulation with a=0.96, and the blue circle line is the 2-flux model simulation with a=0.5. ....................................................................... 29 Figure 4. Profiles of upward TB simulated by the 1-flux and 2-flux models using the same absorption and scattering coefficients for different snow grain size (a), with the equivalent q of 1-flux model calculated from 2-flux model simulations in (b). In (a), q xiv was 0.96 in the 1-flux model simulation (red), where a was calculated by the Improved Born Approximation in the 2-flux model simulation (blue) at 37 GHz. Dots, circles and squares are for small (pec=0.05 mm), middle (pec=0.18 mm) and large (pec=0.32 mm) size of snow particles, respectively. ......................................................................................... 30 Figure 5. Upward TB at the snow surface simulated by the 1-flux and 2-flux models using the same absorption and scattering coefficients by Improved Born Approximation at 37 GHz: (a) TB versus snow slant thickness for medium snow grain size (pec=0.18 mm), (b) TB versus exponential correlation length (pec) for 50 cm snow slant thickness................ 32 Figure 6. Scattering coefficient used by the HUT model and MEMLS as a function of pec at 37 GHz. ......................................................................................................................... 33 Figure 7. Backward-scattering coefficient of MEMLS and 4% of the scattering coefficient used by the HUT model as a function of pec at 37 GHz. ................................ 35 Figure 8. Upward TB at the snow surface at 37 GHz simulated by the 1-flux and 2-flux models using the different scattering coefficients of the HUT model (red lines) and MEMLS, with 6-flux theory included (green lines) and not included (blue lines): (a) TB versus snow slant thickness for medium snow grain size (pec=0.18 mm, correspondingly Dmax =1 mm), (b) TB versus exponential correlation length for 50 cm snow thickness. ... 37 Figure 9. The simulated TB of MEMLS and HUT compared to TB observations, at 37 GHz in vertical polarization, where the black dash line is the 1:1 line. The R2 is the coefficient of determinant between the simulation and the observation (i.e., the R2 of the 1:1 line), and the R2reg is the coefficient of determinant for the fitted lines between the observed and the simulated TB. ......................................................................................... 45 xv Figure 10. Bias of MEMLS and HUT simulations compared to TB observed by radiometers as a function of snow depth at 37 GHz in vertical polarization. The black dash line is where dTB equals zero. .................................................................................. 47 Figure 11. Bias of MEMLS and HUT simulations compared to TB observed by radiometers as a function of snow grain size at 37 GHz in vertical polarization. The black dash line is where dTB equals zero. ................................................................................... 48 Figure 12. Bias of MEMLS and HUT simulations compared to TB observed by radiometers as a function of snow depth at 19 GHz in vertical polarization. The black dash line is where dTB equals zero. .................................................................................. 50 Figure 13. Bias of MEMLS and HUT simulations compared to TB observed by radiometers as a function of snow grain size at 19 GHz in vertical polarization. The black dash line is where dTB equals zero. .................................................................................. 51 Figure 14. Observed SnowScat VV-pol. backscattering coefficients (lines) at Sodankla compared to MEMLS3&a simulations (crossed) at four Intensive Observation Periods (IOP). The blue, red and green colors are for 10.2, 13.3 and 16.7 GHz, respectively...... 59 Figure 15. Observed SnowScat VH-pol. backscattering coefficients at Sodankla compared to MEMLS3&a simuations at four IOPs. ......................................................... 61 Figure 16. Measured backscattering coefficient, snow depth & soil frost depth, soil moisture content and snow grain size (Dmax) (from top to bottom) at Sodankyla during IOP1. The gray vertical lines mark important events from snowcover, soil temperature and frost depth measurements; The red dash lines indicate a change of features for the backscattering coefficient measurement. .......................................................................... 65 xvi Figure 17. Measured backscattering coefficient, snow depth & soil frost depth, soil moisture content and snow grain size (Dmax) (from top to bottom) at Sodankyla during IOP2. ................................................................................................................................. 66 Figure 18. Measured backscattering coefficient, snow depth & soil frost depth, soil moisture content and snow grain size (Dmax) (from top to bottom) at Sodankyla during IOP3. ................................................................................................................................. 67 Figure 19. Measured backscattering coefficient, snow depth & soil frost depth, soil moisture content and snow grain size (Dmax) (from top to bottom) at Sodankyla during IOP4. ................................................................................................................................. 68 Figure 20. A schematic view of the flow chart of the MCMC algorithm at the ith iteration. The iteration i starts from the up left corner, where a candidate combination of snow/soil parameters, θi,candidate is produced using jump functions from θi-1. θ includes multi-layer snow thickness, temperature, density, exponential correlation length, and soil temperature, soil volumetric moisture and soil surface roughness. Then, the posterior probability p(θi,candidate |y) and p(θi-1 |y) were calculated and compared using equation (3.3). The algorithm will determine final estimate at i-th iteration (θi) according to r with relevant to a probabilistic value, rc. .................................................................................................... 76 Figure 21. The in-situ measurements (a) and the priors (b)-(c) for a snowpit at Sodankyla observed on March 2nd, 2010 as an example. The local priors used the average values of the measured snow properties at Sodankyla in March, see (b). The generic priors used the VIC-model SWE predictions and the Sturm’s snow classes; the generic prior of pec is 0.18±0.18 mm, see (c). The priors for layer thickness were converted from SWE and xvii density priors. On the up-right corner, it shows the geographic location of the sites used in this paper. ...................................................................................................................... 83 Figure 22. The performance of the BASE-PM (Bayesian Algorithm for SWE EstimationPassive Microwave) for SWE retrieval using real TB measurements: (a) the BASE-PM SWE compared to the observed SWE at different sites and months, (b) the BASE-PM SWE uncertainty compared to the absolute SWE error. BASE-PM used the generic priors for both SWE and the snow stratigraphy feature. ............................................................. 86 Figure 23. The MCMC-estimated SWE compared to the observed SWE using observed TB from radiometers, when (a) local SWE prior but generic stratigraphy (Strati.) prior were used; (c) generic SWE prior but local stratigraphy prior were used; (e) local SWE prior and local stratigraphy prior were used. In (b), (d) and (f), it shows the MCMC SWE uncertainty versus the absolute SWE error, corresponding to (a), (c) and (e), respectively. ........................................................................................................................................... 88 Figure 24. The MCMC-estimated SWE compared to the observed SWE using synthetic TB, when (a) generic SWE prior and generic stratigraphy (Strati.) prior were used; (b) local SWE prior but generic stratigraphy prior were used; (c) generic SWE prior but local stratigraphy prior were used; (d) local SWE prior and local stratigraphy prior were used. ........................................................................................................................................... 90 Figure 25. The BASE-PM Markov chains of the snowpit observed at Sodankylä on March 2nd, 2010 for (a) snow layer thickness, (b) snow density, (c) pec, and (d) temperature, when generic SWE and stratigraphy priors were used. BASE-PM chose 2 xviii layers, where “BotLyr” is for the chain of the bottom layer, and “SurfLyr” is for the chain of the surface layer. ........................................................................................................... 93 Figure 26. Simulated TB from the estimated snow and soil variables in Figure 19, where the crosses are the measured TB at four frequencies. ........................................................ 94 Figure 27. Histogram of the SWE for the example Sodankylä snowpit, where the blue dash line indicates the BASE-PM SWE, and the back dash line indicates the measured SWE. ................................................................................................................................. 94 Figure 28. Histogram of the chains in Figure 25 after the burn-in period for: (a) layer thickness, (b) density (b), (c) pec, and (d) 274 K minus snow temperature in K. In (d), the variable estimated related to temperature is 274 K minus the temperature in K, because the assumption of a log-normal distribution requires positive values. ............................. 96 Figure 29. TB at 6.925 to 89 GHz versus SWE for all the 48 snowpits, where in (a) it shows the measured TB, in (b) it shows the synthetic TB................................................ 100 Figure 30. Simulated TB at 10.65 to 89 GHz using the measured snow density, temperature and pec for each snowpit but the snow depth was scaled to match the SWE from 20 mm to 400 mm in X-axis. Each line corresponds to the TB-SWE relationship for one snowpit at one frequency: (a) for the snowpits with SWE between 40 and 50 mm; (b) for the snowpits with SWE between 50 and 170 mm; (c) for the snowpits with SWE larger than 170 mm (i.e. the CLPX snowpits). ............................................................... 101 Figure 31. The SWE using 1-layer snow parameter estimation when the same priors as BASE-PM were used: (a) the MCMC SWE versus the true SWE; (b) the estimated SWE uncertainty versus the absolute SWE error. .................................................................... 103 xix Figure 32. Simulated TB at 10.65 to 89 GHz using the 1-layer (first row) or 2-layer (second row) snow density, temperature and pec for each snowpit, where the snow depth was scaled to match the SWE from 20 mm to 400 mm in X-axis. Each line corresponds to the TB-SWE relationship for one snowpit at one frequency. From left to right: the first column is for the snowpits with SWE between 40 and 50 mm; the second column is for the snowpits with SWE between 50 and 170 mm; the third column is for the snowpits with SWE larger than 170 mm (i.e. the CLPX snowpits). .............................................. 104 Figure 33. Simulated upward propagated TB at 89 GHz, vertical polarization for the example snowpit. From left to right, 1-layer, 2-layer and the originally-measured 6-layer snow profiles were utilized. The originally-measured snow depth was used for the black line; for blue dash line and red line, the snow depth was scaled by 2/3 and 1/3, respectively, with relative thickness of layers unchanged. ............................................. 105 Figure 34. Markov chains of layer thickness (dz), exponential correlation length (pec) and snow density for the bottom (Lyr1) and the surface (Lyr2) snow layers (left), and their posterior histogram compared to the prior histogram (right).......................................... 111 Figure 35. Backscattering coefficient calculated from the snow parameters in Markov chains in Figure 34. ......................................................................................................... 112 Figure 36. Sensitivity of backscattering coefficient to the bottom- and surface-layer thickness, pec and density for the example retrieval results in Figure 34. ...................... 114 Figure 37. MCMC-retrieved snow depth, snow density (bulk values) and SWE (from top to bottom) using 3-frequency SnowScat observations at VV pol. based on 2-layer snow estimates.......................................................................................................................... 116 xx Figure 38. Retrieved snow profile for point 1 in Figure 34. ........................................... 117 Figure 39. Retrieved snow profile for point 29 in Figure 34. ......................................... 117 Figure 40. Backscattering coefficient simulated from the MCMC estimates (top) and their errors compared to the SnowScat measurements (bottom)............................................. 118 Figure 41. MCMC-retrieved snow depth, snow density (bulk values) and SWE (from top to bottom) using 3-frequency SnowScat observations at VV pol. based on 1-layer snow estimates.......................................................................................................................... 119 Figure 42. MCMC-retrieved SWE using synthetic-generated backscattering coefficients based on 2-layer snow estimates. .................................................................................... 121 Figure 43. Comparison of retrieved and true bottom- and surface-layer pec and thickness based on 2-layer snow estimates. .................................................................................... 122 Figure 44. MCMC-retrieved SWE using synthetic-generated backscattering coefficients based on 1-layer snow estimates. .................................................................................... 123 Figure 45. Retrieved 1- (top) and 2-layer (bottom) snow profiles for point 9 based on the synthetic-generated backscattering coefficients in Figure 42 and 44. ............................ 124 Figure 46. Retrieved 1- (top) and 2-layer (bottom) snow profiles for point 22 in Figure 42 and 44. ............................................................................................................................. 125 Figure 47. Retrieved 1- (top) and 2-layer (bottom) snow profiles for point 30 in Figure 42 and 44. ............................................................................................................................. 126 Figure 48. Sensitivity of backscattering coefficient to the bottom- and surface-layer thickness, pec and density for point 1 in Figure 37. ........................................................ 128 xxi Chapter 1: Introduction Every year, seasonal snow covers about 50% of the land surface (~46 million km2) in the Northern Hemisphere (Brown & Robinson, 2011). The snow cover reflects up to 80% of the solar radiation, cools the Earth’s surface and impacts the Earth’s meteorological and climatological systems (Flanner et al., 2011). The water resources stored in accumulated snow from winter precipitation provide water supplies for about one-sixth of the world’s population (Barnett et al., 2005). Estimating snow depth, snow mass and the timing of snowmelt are of cardinal importance for hydraulic and hydrological applications (Lettenmaier et al., 2015). In situ snow surveys over large spatial areas are laborious and resource-intensive in polar and high mountainous regions due to extreme environmental conditions. Utilizing airborne and satellite sensors to measure the snow-emitted or reflected electromagnetic waves at appropriate bands with well-designed sensor geometry can remotely collect the information contained in snow via remote sensing techniques. Snow water equivalent (SWE) is the equivalent depth of liquid water contained in snow when snow completely melts (Takala et al., 2011). Electromagnetic waves at microwave bands have wavelengths that are comparable to the dimension of the snow crystals. Electromagnetic radiation at microwave frequencies penetrates through the snow surface, and interacts with the snow crystals via absorption and volume scattering processes. 1 Snow volume scattering strongly reduces the microwave emission from the underlying soil background, and scatters the incident radar radiation to the backscattering direction (Ulaby et al, 2014). Therefore, both the passive microwave brightness temperature (TB) measured by radiometers and actively microwave backscattering coefficient (σ0) measured by radars or scatterometers can be utilized to identify the extent of volume scattering in snowpacks. Such information can then be interpreted by the microwave snow observation models, and related to the physical snow properties. For hydrological and climatological applications, the required accuracy of SWE is ~30 mm for shallow snow (<300 mm SWE) and 10% for deep snow (IGOS, 2007). There have been publically-available passive microwave SWE products (Tedesco et al., 2004a; Takala et al., 2011; Tedesco, 2015). They utilized the brightness temperature observed at about 18.7 and 36.5 GHz from the space-borne SMMR (Scanning Multichannel Microwave Radiometer), SSM/I(S) (Special Sensor Microwave Imager/Sounder), AMSR-E (Advanced Microwave Scanning Radiometer-Earth Observing System) and AMSR2 (Advanced Microwave Scanning Radiometer 2) sensors. The GlobSnow product has an accuracy of 40-mm rms error in flat regions, but it requires snow depth measurements at nearby synoptic stations to provide prior information (Takala et al., 2011). In the past, there have been scientific plans for satellite radar missions, such as the Cold Regions Hydrology High-resolution Observatory (CoReH2O) mission proposed to the European Space Agency (ESA) (ESA, 2012) and the Snow and Cold Land Processes (SCLP) mission proposed to the National Aeronautics and Space Administration (NASA). These missions proposed satellite-borne X- and Ku-band dual polarization 2 radars; the scientific goal is to improve the spatial resolution of the current passive microwave SWE products, which will especially benefit the hydrological applications for terrain-complex mountainous regions. Large ground-based snow experiments, such as the Cold Land Processes Field Experiments (CLPX), the Nordic Snow Radar (NoSREx) Experiments (Lemmetyinen et al., 2016), the Canadian CoReH2O Snow and Ice Experiment, and the on-going SnowEx field campaign (Rodell, 2016), measured the physical properties and the microwave signals of natural snowpacks. The high-quality legacy measurements from these experiments can be used to explore the feasibility to improve the current SWE retrieval accuracy based on advanced retrieval algorithms. The research goals of my dissertation are: (1) to evaluate the accuracy of the snow radiative transfer (RT) models used to simulate the microwave signals from physical snow properties, compare the model theories and practical performance, and evaluate their feasibility to work as an observation model in the SWE retrieval. (2) to develop a SWE retrieval algorithm based on a Bayesian-based Markov Chain Monte Carlo approach which can iteratively find the best SWE estimate that matches the microwave observations, apply this algorithm for passive microwave TB and snowpit measurements from the ground-based experiments, and study the required prior information to support the retrieval. 3 (3) to adapt the retrieval algorithm to active microwave SWE retrieval, evaluate the algorithm performance via synthetic test, and discuss the factors that influence the algorithm performance. The novel points of this research compared to previous studies are: (1) it aims to build a new snow RT model comparison framework that decomposes the current commonly-used models into different aspects of theories and assumptions, and do model comparisons in parallel according to these differences using a single notation system for all model parameters. (2) it aims to consider the vertically-inhomogenous snow properties of natural snowpacks, estimating snow properties for multiple layers simultaneously, and explore whether this strategy is necessary for passive and active SWE retrieval. (3) it aims to use microwave measurements at more frequencies than existing algorithms in order to provide stronger constraints on the SWE retrieval process to improve the retrieval accuracy. In the following sections, progress of the related studies and the research plans for the goals mentioned above are briefly introduced. Section 1.1: Snow Radiative Transfer Models Snow radiative transfer (RT) models estimate microwave radiation that would be observed by a sensor given a set of natural snow properties, which include snow depth (thickness), snow temperature, snow density, and the snow micro-structure, often characterized by grain size. According to Fierz et al. (2009), snow depth is the vertical 4 snowpack thickness. Snow density is defined as the weight of snow in a unit volume (in kg/m3), and the snow density divided by the pure ice density (917 kg/m3) is called the snow volume fraction. Snow is a highly porous, sintered medium composed of ice and air. The dimensions of snow micro-scale inhomogeneity can be described by different snow micro-structure metrics listed in Matzler (2002). The development and evaluation of snow RT model can improve our understanding of the snow emission and scattering physics (Kontu & Pulliainen, 2010; Matzler, 1987; Derksen et al., 2012). A physically based SWE retrieval algorithm requires good forward RT models with well-characterized simulation accuracy and reasonable sensitivity to snow variables (Durand et al., 2008; Takala et al., 2011). Original snow RT models considered snow as a one-layer homogenous medium (Pulliainen et al, 1999; Tsang et al., 2000; Jiang et al, 2007), but hitherto, multiple-layer models have been developed (Lemmetyinen et al., 2010; Wiesmann &Matzler, 1999; Picard et al., 2013). The commonly-used snow models, such as the Helsinki University of Technology (HUT) snow emission model (Lemmetyinen et al., 2010), Microwave Emission Model of Layered Snowpacks (MEMLS) (Wiesmann &Matzler, 1999), Dense Media Radiative Transfer model based on Quasi-Crystalline Approximation (QCA) Mie scattering of Sticky spheres (DMRT-QMS) (Tsang et al., 2000), the bi-continuous model (Ding et al., 2010; Chang et al., 2014) and Dense Media Radiative Transfer - Multi Layer model (DMRT-ML) (Picard et al., 2013), differ in snow micro-structure assumptions, snow volume scattering theories, and radiation flux angular resolution used in the RT solver. To describe the micro-structure of snow, DMRT models treat snow as aggregates 5 of discrete ice spheres (Tsang et al., 1992; Picard et al., 2013), or clusters of spheres with stickiness between them (Tsang et al., 2000), whereas MEMLS and the bi-continuous model treat snow as a continuous medium. There are three branches of methods to calculate the volume scattering coefficient of snow, in general: the semi-empirical methodology (Pulliainen et al., 1999; Wiesmann & Matzler, 1999), the Monte Carlo simulation (Ding et al., 2010; Chen et al., 2003), and analytical methods (Tsang et al., 1992; Tsang et al., 2000; Wang et al., 2000). Of analytical methods, the strong fluctuation theory (SFT) (Wang et al., 2000) and the Improved Born Approximation (IBA) (Matzler & Wiesmann, 1999) were used for most continuous-medium assumptions; where the QCA theory was used for the ice-particle assumptions (Tsang et al., 1992; Tsang et al., 2002). However, the bi-continuous model, which also considers dry snow as a twocomponent continuous medium of ice and air, used Monte Carlo simulation as an exact solution to Maxwell’s equations (Ding et al., 2010; Chang et al., 2014). To solve the radiative transfer (RT) equation, there are options to split the radiation in the 4π space into two fluxes (Boyarskii & Etkin, 1994), six fluxes (Wiesmann & Matzler, 1999), or multiple fluxes using Gaussian quadrature and Eigen analysis (Tsang et al., 2000; Picard et al., 2013; Wang et al., 2000). The situation becomes more complicated when considering the range of parameters that are used to describe snow-micro-structure. In general, the models based on continuous-medium assumption are more likely to use some statistical parameter, such as the exponential correlation length (pec) for the SFT (Wang et al., 2000) and IBA theories (Matzler & Wiesmann, 1999), and the mean and the standard deviation of wave number in the bi-continuous model (Ding et al., 2010; Chang et al., 6 2014). QCA theory used the diameter of single spheres, as well as a “stickiness” parameter to describe the possibility of spheres sticking together (Tsang et al., 2000). The HUT model uses an effective diameter of snow grains (Pulliainen et al., 1999), which is empirically parameterized from the observed snow grain size (Dmax) defined as the maximum dimension of characteristic snow grains (Matzler, 1997). Due to the errors in the snowpit measurements and the microwave measurements, as well as the different snow features in the datasets used for previous model validations, it is very difficult to evaluate which model is more accurate. Also, because each commonlyused model takes a combination of choices from the model theory aspects mentioned above, it is very difficult the evaluate the source of errors. Early model comparison works only compared the model simulations with the microwave measurements (Tedesco & Kim, 2006a). Follow-on work attempted to address the details in model theories and assumptions. For example, Lowe & Picard (2015) compared the volume scattering coefficient equations of IBA and QCA-CP (Coherent Potential), and proved their similarities in theory and simulations. Matzler (2006) compared the 1st order, 2nd order and high-order iterative methods for RT equation resolution. Roy et al. (2013) studied the differences and the relations in snow micro-structure parameters based on the measurements in natural snowpits; Royer et al. (2017) applied these measurements to HUT, MEMLS, DMRT-QMS and DMRT-ML, and discussed the required scaling factors when they are used as model inputs. In this research, a RT model with moderate calculation cost is considered to be used as the forward observation model in an iterative SWE retrieval system. To be utilized for 7 satellite product in future, it is good to start with simple models. However, we hypothesize that a multi-layer RT equation is mandatory to consider the verticallyinhomogeneous snow properties, and a physically-based or at least semi-empirical volume scattering coefficient is required for simulating microwave signals with good accuracy for a large range of frequencies. Here, the Helsinki University of Technology (HUT) model (Lemmetyinen et al., 2010) and the Microwave Emission Model of Layered Snowpacks (MEMLS) (Wiesmann & Matzler, 1999; Matzler & Wiesmann, 1999) are considered. They have many features in common: for example, they both consider multi-layer snow, were constructed with radiometric measurements of natural snow samples, and used analytical approaches to solve the RT equation which can be stated clearly by comparison with the 2-flux theory. The HUT model has two semi-empirical extinction coefficient parameterizations; MEMLS has a semi-empirical equation and the IBA to calculate the volume scattering coefficient. The differences of the model theories were compared, and the model performance was evaluated using in-situ microwave and snowpit measurements. This comparison is presented in Chapter 2. Section 1.2: Passive SWE Retrieval Algorithms The theory of SWE retrieval from passive microwave brightness temperature (TB) measurements relies on the volume scattering produced by snow crystals (Ulaby and Long, 2014). The volume scattering effect strongly decreases the TB observed the snow surface at about 37 GHz (Rott, 1986), and the number of snow crystals is related to the snow depth or SWE. However, SWE retrieval is not a easy task, not least because the size 8 of the snow crystals also has a strong influence on the magnitude of volume scattering (Armstrong et al., 1993; Kelly et al., 2003). Indeed, according to snow RT models (Wiesmann and Matzler, 1999; Liang et al., 2008; Lemmetyinen et al., 2010; Picard et al., 2013), TB observed at the snow surface is a complicated, non-linear function of a large number of parameters, including snow depth, snow grain size, snow density, snow temperature, and the background emission from the underlying soil determined by soil moisture, temperature and roughness. Snow properties are not vertically homogenous due to several physical processes continuously ongoing in the natural snowpacks (Fierz et al., 2009; Jordan, 1991), and this fact has been confirmed by many ground snowpit measurements (Hardy et al., 2003; Lemmetyinen et a., 2016). Snow stratigraphy has nonnegligible effects on the TB observations (Durand et al., 2011; Pan et al., 2012). Different layer combinations of the snow properties can also produce similar microwave signals, which further complicates the inversion problem. Traditional inversion algorithms are based on empirical relationships between TB and SWE. It was first assumed that TB at high frequencies are sensitive to the amount of scattering in the snow medium, whereas TB at low frequencies provide a background soilemitted TB that is less sensitive to snow. Therefore, the TB difference between low and high frequencies can be related to the snow depth or SWE (Chang et al., 1987). These empirical relationships can be adjusted according to different snow types and the influence of the surrounding vegetation. The differences in band choice and coefficient determination result in a series of static and dynamic algorithms, as summarized in Clifford (2010), Frei et al. (2012) and Jiang et al. (2014). More sophisticated SWE 9 estimation algorithms have since been developed, which rely on the snow RT models to find the best estimate of snow properties that match the observed TB, or to constrain the simulations of the land surface model in a probabilistic-based dynamic system. Details of the data assimilation methods can be found in Andreadis and Lettenmaier (2012). Of all the previous SWE estimation algorithms, the physically-based retrieval algorithms using only microwave measurements can be classified into search-based methods and signal-decomposition methods. In search-based methods, TB signals are simulated using random combinations of snow parameters, and the best combination, which gives a simulated TB signal closest to the observations, is chosen as the retrieval result. The accuracy and the efficiency of this procedure are determined by the algorithm utilized. Some examples are the simple searching methods (Pulliainen et al., 1999; Roy et al., 2004; Butt and Kelly, 2008; Butt, 2009), simple searching methods with prior information (Pulliainen, 2006; Takala et al., 2011), Artificial Neural Network (ANN) methods (Davis et al., 1993; Tedesco et al., 2004b; Gan et al., 2009), and Genetic Algorithm (GA) (Tedesco & Kim, 2006b). In these methods, the one-layer HUT model (Pulliainen et al., 1999) and the one-layer DMRT-QCA model with coherent potential (CP) (Tsang & Kong, 2001) were used as the observation model (Pulliainen, 2006; Tedesco & Kim, 2006b; Takala et al., 2011). Prior information has been used in SWE retrieval. For example, Takala et al. (2011) estimates an effective snow grain size at station locations where the snow depth is known and interpolates it to nearby pixels as the grain size prior. The signal-decomposition methods utilize the empirical relationships found in physically-based snow RT model simulation datasets. Jiang et al. (2007) used 10 two-frequency, dual-polarization TB to separate the snow volume scattering contributed emission and the soil contributed emission. Later, SWE was empirically fitted as a function of the snow-contributed emission. Only one snow layer was considered in Jiang et al. (2007). In this research, the goal is to build a physically-based SWE retrieval algorithm allowing for multi-layer snowpack and prior information on the snow properties. In general, it can be classified as a search-based method with priors. In this algorithm, the Markov Chain Monte Carlo (MCMC) method (Gelman et al., 1995) was applied to choose snow parameters that give the closest simulated TB as the observed values, while taking the prior information into account in a Bayesian sense. It is called the Bayesian Algorithm for SWE Estimation with Passive Microwave measurements (BASE-PM). BASE-PM applies a similar SWE retrieval method as that introduced in Durand and Liu (2012), but it uses actual TB observations from ground-based radiometers, rather than synthetically generated TB estimates. As mentioned in Durand and Liu (2012), prior information is required to prevent the estimated snow variables from going beyond their reasonable range in natural conditions. Therefore, another focus of this research is to how to produce the appropriate prior for practical SWE retrieval. The influence of different priors will be discussed. This part of the research content will be shown in Chapter 3. Section 1.3: Active SWE Retrieval Algorithm The snow volume scattering effect increases the backscattering coefficient measured by active sensors at about 15 to 30 GHz, which can be used to estimate SWE on the ground 11 (Rott, 1986). Backscattering coefficient at 13.9 GHz measured by POLSCAT (POLarimetric SCATterometer) during the CLPX experiment in 2006-2008 increased about 0.15-0.5 dB for every 10 mm increase of SWE (Yueh et al., 2009). The frequency is lower than the Ka-band used for passive microwave SWE retrieval, because the radar pulses have much stronger power than the naturally-emitted microwave radiation, and because radar has stronger interaction with the snow crystals. Both search-based algorithm and signal decomposition algorithm have been used for the active microwave SWE retrieval. Rott et al. (2012) used the semi-empirical Radiative Transfer (sRT) model and a cost function to iteratively find the estimates of SWE and an effective grain size that match the backscattering coefficient observations. Several examples of the signal decomposition algorithm applied for radar SWE retrieval can be found in Shi & Dozier (2000), Du et al. (2010a), Xiong et al. (2014) and Cui et al. (2016). The idea of this algorithm is also to separate the soil-contributed backscattering signals from the total backscattering coefficient observed at the snow surface. Empirical relationships were built from1-layer DMRT-QMS or bi-continuous model generated datasets. In Shi & Dozier (2000) when the C- and X-band copolariaztion (VV & HH) radars were utilized, the soil-contributed backscattering signals were derived from the relations with the corresponding component from the L-band SAR measurements. Du et al. (2010a) applied the relations of both soil- and snow-contributed backscattering signals at two angles to estimate SWE from multi-angle X-band SAR measurements. The algorithms of Xiong et al. (2014) and Cui et al. (2016) were developed for X- and Kuband dual-polarization (VV & VH) SAR. The snow-contributed backscattering 12 coefficient was fitted as a function of VV-pol. backscattering coefficient and the difference of VV- and VH-pol. backscattering coefficients. In general, the signal decomposition algorithms highlight the importance of the soil background, because its backscattering contribution represents a large portion of the observations. In this research, the active SWE retrieval algorithm is adapted from the passive microwave BASE-PM algorithm. It is a follow-on research based the same Markov Chain Monte Carlo method, but the TB observation model is changed to the newly-built Microwave Emission Model of Layered Snowpacks 3&active (MEMLS3&a) (Proksch et al., 2015). The SWE retrieval performance will be tested for both real and syntheticallygenerated backscattering coefficients. This part of the research content will be shown in Chapter 4. 13 Chapter 2: Snow Radiative Transfer Model Comparison In this chapter, the theories and the performance of the Helsinki University of Technology (HUT) snow emission model and the Microwave Emission Model of Layered Snowpacks (MEMLS) are compared. The goal is to highlight model similarities and flag out their differences both in theory and in practice. The approach is firstly to present a concise, side-by-side derivation of model equations using a single notation system. Secondly, simplified, side-by-side experiments were conducted to explore the ways that the differences in theory influence TB simulations. Thirdly, the model performance was evaluated for a large number of snowpits from Steamboat Springs, and Fraser, USA; Churchill, Canada; and Sodankylä, Finland, and the results were explained in terms of the theoretical differences found. Section 2.1: Differences in Model Theory In this section, the major theoretical differences between HUT and MEMLS are presented. The first sub-section shows the derivation of the RT equation from the same general form, and describes where the two models diverge with regard to the 1-flux assumption invoked by the HUT model. The second sub-section describes the differences in the volume scattering coefficient calculation. The third sub-section describes how the trapped radiation is considered in MEMLS by a 6-flux formulation. 14 One-Flux Versus Two-Flux Theories Based on general RT theory, as microwave propagates a homogenous and isotropic snow media, the change of TB in the propagation direction of θ (zenith angle) and φ (azimuth angle), at any vertical location (z) can be written as (Ulaby et al., 2014): ∂TB ( z, θ , ϕ ) 1 = −κ eTB ( z, θ , ϕ ) + κ aT + ∂z ⋅ secθ 4π ∫∫π κ bi s (θ , ϕ , θ ' , ϕ ' )TB ( z, θ ' , ϕ ' ) sin θ ' dθ ' dϕ ' (2.1) 4 where κe is the extinction coefficient, and the first term on the right-hand side (R.H.S.) is the attenuation of microwave due to absorption and scattering; κa is the absorption coefficient, and T is the physical temperature of snow, where the second term on the R.H.S. is the emitted radiation from snow as a heat source; κsbi(θ,φ,θ’,φ’) is the bistatic scattering coefficient, and its averaged value over 4π space is the scattering coefficient, κs: 1 π 2 π bi κs = ∫ ∫ κ (θ ,φ ,θ ',φ ') sin θ 'dθ 'dφ ' 4π θ '=0 φ '=0 s (2.2) where, the sum of κs and κa is κe. The third term on the R.H.S. is the sum of scattered microwaves from other directions into the direction of (θ, φ). If the snow medium is assumed plane-parallel, equation (2.1) can be simplified by reducing TB into two directions as: the “upward” direction (θ0, φ0), and the “downward” direction (θ0+π, φ0+π). Using 2 fluxes, equation (2.1) can be rewritten as: ∂Tup ( z ' ) ∂z ' − = −κ e Tup ( z ' ) + κ a T + κ sforward Tup ( z ' ) + κ sbackward Tdn ( z ' ) ∂Tdn ( z ' ) = −κ eTdn ( z ' ) + κ aT + κ sforward Tdn ( z ' ) + κ sbackward Tup ( z ' ) ∂z ' 15 (2.3) (2.4) where z’=z·secθ0, and θ0 is the observing zenith angle; Tup(z’) is the upward- propagated TB; Tdn(z’) is the downward-propagated TB; κsforward and κsbackward are the 2-flux forward scattering coefficient and backward scattering coefficient, respectively, and the sum of κsforward and κsbackward is κs. If the medium is locally homogenous and there is no refraction due to inhomogeneity, κsforward and κsbackward can be defined as (Matzler & Wiesmann, 1999): π /2 κ sforward = ∫ sin θ dθ θ =0 κ sbackward = π /2 1 π /2 2 π bi ∫ ∫ κ (θ ,φ ,θ ',φ ')sin θ 'dθ 'd φ ' 4π θ '=0 φ '=0 s 1 π 2π ∫ sin θ dθ 4π ∫ ∫ κ bi s (θ , φ ,θ ', φ ')sin θ 'dθ 'd φ ' (2.6) θ '=π /2 φ '=0 θ =0 (2.5) Using κe=κa+κs and κs= κsforward + κsbackward, (2.3) and (2.4) can be rearranged to: ∂Tup (z ') = −κ a Tup (z ') − T − κ sbackward Tup (z ') − Tdn (z ') ) ∂Tdn (z ') = −κ a Tdn (z ') − T − κ sbackward Tdn (z ') − Tup (z ') ∂z ' ) ∂z ' − ) ( ( ) ( ( (2.7) (2.8) These are the differential RT equations used in MEMLS (Wiesmann & Matzler, 1999a). Because the upward TB and downward TB are interlinked, it is called the “2-flux” model in this paper. The forward scattering coefficient was eliminated algebraically; therefore, it actually only uses the backward scattering coefficient. The HUT model assumes that the scattered intensity is mostly concentrated in the forward direction, and introduces an empirical parameter, q, to simplify equation (2) into: ∂Tup (z ') ∂z ' = κ aT + (qκ s − κ e )Tup (z ') (2.9) 16 − ∂Tdn (z ') = κ aT + (qκ s − κ e )Tdn (z ') ∂z ' (2.10) The key point of this simplification is the sum of the scattered intensities from all directions can be empirically expressed by the scattered intensity in the forward direction scaled by a ratio of q. In this case, the change of upward TB depends only on upward TB, thus we call it “1-flux” model. It does not mean in the HUT model, the upward and downward radiations are entirely independent. The reflection at layer interfaces reflects the upward fluxes into downward fluxes, or downward fluxes into upward fluxes. The transmission at layer interfaces links the fluxes between different adjacent layers. However, inside a homogenous medium, it is really a “1-flux” model. By assuming the upward TB of the 1-flux and the 2-flux models to be equal between equation (2.3) and (2.9), it yields: q( z ' ) = κ sforward Tup ( z ' ) + κ sbackwardTdn ( z ' ) T ( z' ) = a + (1 − a) ⋅ dn κ sTup ( z ' ) Tup ( z ' ) (2.11) where, κ sforward a= κs (2.12) q here is the forward scattered ratio of intensities, where a is forward-scattered ratio of scattering coefficients. They are not the same unless the scattered intensities in the backward direction (κsbackwardTdn(z’)) equals zero. Using the experimental data in Matzler (1987) and Hallikainen et al.(1987), q in the HUT model was estimated to be 0.96 for all snow types, all frequencies, and at any vertical location inside the snow cover. 17 z 0 d,T -d θ T up (0 − ) T dn (0 − ) T up (− d + ) T dn (− d + ) Figure 1. A homogenous snowpack of thickness, d, where z=0 is at the snow surface, and z=-d is at the bottom of snow. For homogenous snow layer (1-layer) as shown in Figure 1, the general solution for the second-order differential equations (2.7)-(2.8) in the 2-flux model is: Tdn (z ') = T + A0 exp(γ z ') + B0 exp(−γ z ') Tup (z ') = T + γ 0 A0 exp(γ z ') + (2.13) 1 B exp(−γ z ') γ0 0 (2.14) where, κ sbackward (1− a)κ s = κ a + κ sbackward + γ κ a + (1− a)κ s + γ (2.15) γ = κ a2 + 2κ aκ sbackward = κ a2 + 2(1− a)κ aκ s (2.16) γ0 = The constants, A0 and B0, need to be determined by using boundary conditions. To avoid the complexity of considering reflections at layer interfaces, here the downward TB just below the snow surface, Tdn(0-), and the upward TB just above the bottom of the layer, Tup(-d’+), are used. Then, A0 and B0 could be solved as shown in Appendix A. 18 Substituting A0 and B0 into equation (2.13) and (2.14) yields: Tup (z ') = T + 1 " T (−d '+ ) − T η (−z ') + Tdn (0− ) − T ζ (z '+ d ')$% η (d ') # up (2.17) Tdn (z ') = T + 1 " T (0− ) − T η (z '+ d ') + Tup (−d '+ ) − T ζ (−z ')$% η (d ') # dn (2.18) ( ( ) ) ( ( ) ) where, d’ = d·secθ0; d is the snow depth, and secθ0 is the observing angle. η and ζ are two intermediate functions defined as: η (x) = 1 exp(γ x) − γ o exp(−γ x) γ0 (2.19) ζ (x) = exp(γ x) − exp(−γ x) (2.20) where, the independent variable x represents certain amount of snow slant thickness inside the snowpack. In the 1-flux model, only the lower boundary condition Tup(-d’+) is needed to solve for the upward TB. Likewise, only the upper boundary condition Tdn(0-) is needed to solve for the downward TB. The solutions to equation (2.9) and (2.10), respectively, are: Tup ( z ' ) = κ aT (1 − exp[−(κ e − qκ s )(z'+d ' )]) + Tup (−d '+ ) ⋅ exp[−(κ e − qκ s )(z'+d ' )] (2.21) κ e − qκ s Tdn ( z ' ) = κ aT (1 − exp[(κ e − qκ s ) z ' ]) + Tdn (0 − ) exp[(κ e − qκ s ) z ' ] κ e − qκ s (2.22) Comparing (2.17) and (2.21), it shows the upward TB is still a function of the downwardpropagated boundary condition in the 2-flux model, but not in the 1-flux model. 19 The difference between the 1-flux and 2-flux models can be better understood by considering a snowpack with infinite thickness. When d’ approaches infinity, the upward TB just below the snow surface will approach: Tup，1− flux (0 − , d ' → +∞) → κa κ e − qκ s (2.23) T Tup，2− flux (0 − , d ' → ∞) → T (1 − γ 0 ) + Tdn (0 − )γ 0 (2.24) This simplified case shows that, in equation (2.24), γ0 could be considered as the reflectivity and 1-γ0 as the emissivity of an infinite snow medium, which is more or less named the same as in Wiesmann & Matzler (1999a). However, although q was defined as the forward-scattered ratio of intensities with a clear physical meaning, equation (2.23) has the coefficients in front of the temperatures on the right-hand side, which do not sum up to 1. It thus violates the Kirchhoff’s Law of thermal radiation. To solve the RT equation for multiple-layer snow, it is more convenient to make use of the four radiation quantities as show in Figure 2 (Wiesmann & Matzler, 1999a), where Aj and Bj are downward TB and upward TB just above the lower boundary of layer j, respectively; and, Cj and Dj are the downward TB and upward TB just below the upper boundary of layer j, respectively. Appendix B shows how to solve the 1-flux and 2-flux models for multiple-layer snow. In general, they link the upward and downward TB between adjacent snow layers by the Fresnel equations (assuming smooth snow layer interfaces), and use the RT equation for homogenous snow to link the fluxes inside each layer. In the second process, the difference between the 1-flux and 2-flux models is introduced. By reproducing the solutions that are exactly the same as in Lemmetyinen et 20 al. (2010) for the multiple-layer HUT model, and in Wiesmann & Matzler (1999a) for MEMLS, respectively, we managed to prove that the very different apparent forms of the multiple-layer RT equations are due only to the difference between the 1-flux and 2-flux assumptions. In other words, the differences can be tracked down to the basic differential RT equations shown in equation (2.7)-(2.8) and (2.9)-(2.10). It also shows the different terms of radiations in the HUT model and MEMLS are rigorously given and physically solved from the RT equations. z zn = 0 zj z j-1 An+1 = Ts k y . . . dj , T j . . . z0 =-d Cn Dn A j+1 B j+1 Cj Aj C j-1 Dj Bj D j-1 A1 B1 sn = s as sj s j-1 D 0 = Tg r o u n d s0 = sgs Figure 2. Four TB quantities, Aj, Bj, Cj and Dj, inside the j-th snow layer defined in Wiesmann & Matzler (1999a). The z-axis is along the vertically upward direction. zj-1 and zj indicate the lower and upper boundary of the j-th snow layer. In the figure, sgs is the ground-snow boundary reflectivity; sas is the air-snow boundary reflectivity; Tground is the upward ground TB; and, Tsky is the downward sky TB. Scattering Coefficient Both MEMLS and HUT require extinction coefficients. These quantities are conceptually identical between the two models. 21 Both the HUT model and MEMLS calculate κa from the snow permittivity: κ a [1/m] = 4π ⋅ img ( ε) s λ0 (2.25) where, εs is the snow permittivity, with its real part determined by snow density, and its imaginary part determined by snow density and temperature (Wiesmann & Matzler, 1999); img() means to take the imaginary part of the term inside parenthesis; λ0 is the vacuum wavelength. A preliminary difference between the HUT and MEMLS κs is their approach to calculate volume scattering with different metric to describe snow microstructure. The HUT model parameterizes κs using snow grain size, but MEMLS uses exponential correlation length (pec). MEMLS assumes the snow is a two-phase continuous media. In Matzler (2002), pec is defined as the parameter to fit the spatial auto-correlation function of the snow medium as A(x)=exp(-x/pec). HUT model assumes snow is composed by irregular discrete ice particles. For the HUT model, Hallikainen et al. (1987) measured the power transmission loss of natural snow samples, and derived an empirical equation for κe as: ⎧ ⎪ 0.0018 f 2.8 D02.0 for 18-60 GHz κ e [1 / dB] = ⎨ 300D01.9 for 90 GHz ⎪⎩ (2.26) where, D0 is the snow grain size (mm) and f is frequency (GHz). The equation is thereby called the “Hallikainen equation”. Roy et al. (2004) determined an equation for κe by empirically fitting one-layer HUT model with airborne-observed TB at 18 and 37 GHz, as: 22 κ e [1 / dB]=2 ( f 4 D06 ) 0.2 (2.27) This equation is called the “Roy equation”. The κe in 1/dB in these two equations can be compared directly to the MEMLS values given below by converting to units of 1/m by dividing 10 log10 e = 4.3429 (Rees, 2001). When the HUT model was applied for the TB simulations, all measured snow grain size (Dmax), as “the average of the maximum extension (diameter) of the prevailing or characteristic grains”, were adjusted using an effective value, Deff, as (Kontu & Pulliainen, 2010; Lemmetyinen et al., 2010): D0 [mm] = Deff [mm] = 1.5[mm]⋅ (1− exp(−1.5⋅ Dmax [mm])) (2.28) The method aims to reduce the variations in the large values of the observed grain size, and is taken as one part of the data processing procedures used both for the Hallikainen and Roy equations. For MEMLS, Wiesmann et al. (1998) measured TB of snow slabs (pec from 0.03 to 0.32 mm, density from 100 to 400 k/gm3) over an absorber and a metal plate at 11, 21, 35, 48 and 94 GHz, and solved the scattering coefficient by inverting a one-layer MEMLS. The retrieved κs was fitted as (Wiesmann & Matzler, 1999a), 2.5 2.5 ' ! p $ * ! $ pec f ec ) , κ s [1/m] = 3.16 + 295# & ⋅# & ) 1 mm 1 mm % , " 50 GHz % " ( + (2.29) where pec is the exponential correlation length (mm) and f is the frequency (GHz). It should be noted that a was considered to be 0.5 in this empirical method, because solely 23 by inverting MEMLS formulated with a backward-scattering coefficient, there is no way to acquire any information about the scattering in forward direction. Mätzler (1998) and Mätzler & Wiesmann (1999) derived and tested the improved Born approximation (IBA), where the bistatic volume-scattering coefficient is calculated as, κ sbi (θ , φ ,θ ', φ ') = v(1− v)(εi −1) 2 K 2 Ik 4 sin 2 χ (2.30) where, (θ’, φ’) is the incident angle; (θ, φ) is the scattered angle; v is the ice volume fraction; εi is the ice permittivity; K2 is the ratio of the mean-squared electric fields inside and outside of ice particles; I is defined as an integral of the spatial autocorrelation function, thus it is a function of pec; and, χ is the angle between the direction of the electric field of the incident wave and the propagating direction of the scattered wave. Trapped Radiation In solving RT equation in 3-D space, MEMLS considers the internal reflection in snow medium due to the higher permittivity of snow compared to air. The internal reflection traps the radiation along a subset of directions inside the snow cover, and MEMLS divides the trapped radiation into four horizontal fluxes (Wiesmann & Matzler, 1999a). The gradient of these horizontal fluxes is zero, because snow is assumed plane-parallel in MEMLS (Matzler, 2004b). Therefore, together with the upward and downward fluxes, the radiation in 3-D space is divided into six fluxes in total, and accordingly six RT equations are constructed. These six RT equations can be reduced to a pair of 2-flux RT equations, but with absorption and backward-scattering coefficients changed to: 24 ⎛ κ scross ' ⎞ κ a,6-flux = κ a ⎜1+ 4 ⎟ κ a + 2κ scross ' ⎠ ⎝ backward κ s,6-flux = κ sbackward '+ 4 (2.31) κ scross '2 κ a + 2κ scross ' (2.32) where, κsbackward’ and κscross’ (γb and γc in (Wiesmann & Matzler, 1999a)) are 6-flux backward scattering coefficient and cross scattering coefficient, respectively. The difference between κsbackward and κsbackward’ is that, κsbackward is the integral from bistatic scattering coefficient in the backward-scattered half space, whereas the integral of κsbackward’excludes the bistatic scattering coefficient beyond the critical angle. Section 2.2: Experiments to Test Three Model Theoretical Differences As mentioned in the previous section, there are three major theoretical differences between MEMLS and the HUT model: the 1-flux versus 2-flux theory, the volume scattering coefficient theory, and the trapped radiation. In addition to these three factors, The HUT model and MEMLS also have other differences: MEMLS considers the polarization mixing and the coherency of the reflected radiation by the top and the bottom of a very thin snow layer. However, the effect of polarization mixing is minor for seasonal snow cover (less than 1 K). Coherent effect appears for only part of the snowpits. Therefore, these two points are not included in analysis of model theoretical differences. In this section, numerical experiments are designed to test the effects of three model theoretical differences. 25 Experiment Design 1) One-Flux Versus Two-Flux Theories The goal of the first set of experiments is to evaluate the effect of the 1-flux assumption used in HUT vs. the more general 2-flux theory used in MEMLS. Rather than running the full model codes, the 1-flux and 2-flux models were used to calculate the upward TB just below the snow surface, Tup(0-), using the same absorption and scattering coefficients. Some representative cases were chosen instead of covering an exhaustive range of snow properties. An absorption coefficient of 0.226 1/m was used for both models, which was calculated from a snow density of 200 kg/3 and a snow temperature of -5 °C at 37 GHz. Three sizes of snow particles were considered for the scattering coefficients, as a pec (exponential correlation length) of 0.05 mm, 0.18 mm and 0.32 mm; and the values at 37 GHz were calculated by the physically-based Improved Born Approximation (IBA), as listed in Table 1. Boundary conditions of Tdn(0-) as 2.7 K and Tup(-d+) as 265 K were used throughout the experiments in Section 2.2. Tup(-d+) is much larger than Tdn(0-) because for real snowpacks, Tdn(0-) represents the downwelling radiation from the sky radiation, and Tdn(0-) represents the upwelling radiation from the underlying soil beneath the snow. Size of Snow Particles pec (mm) Scattering coefficient (1/m), κs Small Middle Large 0.05 0.18 0.32 0.0210 0.897 4.202 Backwardscattering coefficient (1/m), κsbackward 0.0105 0.435 1.920 Forward-tototal ratio of scattering coefficients, a 0.501 0.514 0.543 Asymmetry factor, g 6.73*10-5 0.035 0.475 Table 1. The scattering coefficients calculated by the Improved Born Approximation at 37 GHz for three sizes of snow particles 26 2) Scattering Coefficient To compare the scattering coefficients of HUT and MEMLS, a relationship between the snow grain size (Dmax) used by HUT and the exponential correlation length (pec) used by MEMLS is required, and here the method in Durand et al. (2008) was used: ! a + a ln D [mm]± ε , if v > v & D > d # 0 1 max fit 0 max 0 pec [mm] = " p [mm]± ε , otherwise #$ 0 0 (2.33) where, a0 and a1 are two empirical factors, v0 is the threshold for small snow volume fraction, d0 is the threshold for small snow grain size. εfit and ε0 are the fitted error and the uncertainty. The method was built by analyzing the micro-structure of snow samples from Weissflühjoch, Davos, Switzerland in Matzler (2002). In the classification system of Sturm et al. (1995), the snow at Weissflühjoch is alpine snow. The dataset had a wide range of pec (0.034-0.33 mm) and density (93-345 kg m-3) values, and as it was found that there is good empirical relationship between pec and the natural logarithm of Dmax (Durand et al., 2008). For low-density snow, the pec remains a low value as p0 that is unrelated to Dmax. As an example of this point, consider newly-fallen snow, which can have a large geometric dimension, but keeps a low density and low correlation due to its dendritic structure. The data at Weissflühjoch fits a0 as 0.18 mm, a1 as 0.09 mm, and p0 as 0.05 mm; the thresholds are v0 as 0.2 and d0 as 0.125 mm. εfit is 0.027 mm, and ε0 is 0.017 mm. The HUT-model κs was calculated for pec from 0.05 mm to 0.32 mm, and the MEMLS κs was calculated for the corresponding range of Dmax from 0.24 mm to 4.74 mm. Then, the 27 different κs were used in the 1-flux and 2-flux models, respectively, to compare the upward TB at the snow surface, Tup(0-). 3) Trapped Radiation In the last experiment, the effect of 6-flux theory on MEMLS simulations was simulated by comparing the TB at the snow surface before and after applying 6-flux theory. Simulated TB Differences from Model Theories 1) One-Flux Versus Two-Flux Theories Figure 3 compares the TB profile of 1-flux and 2-flux models using the same scattering coefficient (0.897 1/m) for 0.5 and 0.96 a and q values. It should be noted that, this is a sensitivity test, and the fundamental connections between κs and the forward-scattering ratios has not been taken into consideration (i.e., usually a large κs comes with a large a). The value of 0.891 1/m refers to IBA for pec=0.18 mm at 37 GHz. The slant thickness (d’=d secθ0) here is 50 cm. The graphs represent direct solutions to equations (2.17), (2.18), (2.21), and (2.22), and they are not outputs from MEMLS or the HUT model, per se. First, as one could expect, the strongly forward scattering case with a=0.96 gives larger upward TB than the symmetric forward and backward scattering case with a=0.5; so is the situation for 1-flux model, where q=0.96 gives larger TB than q=0.5. It is because one important source of the radiation is the Tup(-d+) beneath the snow, and the forward scattering of upward ground radiation contributes to the upward TB at the snow surface. 28 A comparison between the two scattering cases shows that the 1-flux and 2-flux models give closer TB profiles for strongly forward-scattering case (a=0.96, and q=0.96) than the symmetric forward and backward scattering case (a=0.5, and q=0.5). The former gives a difference of TB at snow surface less than 0.1 K, where the later gives a difference of about 7 K. It is because according to equation (2.11), the difference between the equivalent q of the 1-flux model and a of the 2-flux model is (1-a) multiplied by Tdn(z’)/Tup(z’). Tdn(z’) is smaller in a strongly forward-scattered case. Therefore, the 1flux simplified RT theory is a better approximation to the 2-flux theory in a stronglyforward scattering media. 0 z’ (cm) −10 −20 −30 −40 −50 200 2−flux(a=0.5) 2−flux(a=0.96) 1−flux(q=0.5) 1−flux(q=0.96) 220 240 TBup (K) 260 280 Figure 3. Upward TB profiles of the 1-flux and 2-flux models using the same absorption coefficient and scattering coefficient. The value of scattering coefficient refers to the Improved Born Approximation (IBA) for snow with medium grain size (pec=0.18 mm) at 37 GHz, and the absorption coefficient refers to snow with a density of 200 kg m-3 and a physical temperature of -5 °C at 37 GHz. The red thick line is the 1-flux model simulation with q=0.96. The red thick dash line is the 1-flux model simulation with q=0.5. The blue square line is the 2-flux model simulation with a=0.96, and the blue circle line is the 2-flux model simulation with a=0.5. 29 (a) −10 −10 −20 −20 −30 −40 1−flux(Small) 1−flux(Middle) 1−flux(Large) 2−flux(Small) 2−flux(Middle) 2−flux(Large) Small Middle Large −30 −40 −50 −50 140 160 180 200 220 240 260 280 0.5 TBup (K) (b) 0 z’ (cm) z’ (cm) 0 0.6 0.7 0.8 0.9 q solved from 2−flux model Figure 4. Profiles of upward TB simulated by the 1-flux and 2-flux models using the same absorption and scattering coefficients for different snow grain size (a), with the equivalent q of 1-flux model calculated from 2-flux model simulations in (b). In (a), q was 0.96 in the 1-flux model simulation (red), where a was calculated by the Improved Born Approximation in the 2-flux model simulation (blue) at 37 GHz. Dots, circles and squares are for small (pec=0.05 mm), middle (pec=0.18 mm) and large (pec=0.32 mm) size of snow particles, respectively. Figure 4 shows the performance of the 2-flux and 1-flux models at 37 GHz with identical κs but their default values of a and q. The κs and a used by the 2-flux model refer to IBA for three snow grain sizes described by pec. For the 1-flux model, the same values of κs are used, but q is fixed as 0.96 as in the HUT model. In Figure 4(a), it shows if q is fixed to 0.96, it will produce 1 K, 38 K and 97 K differences in the simulated TB at the snow surface between 1- and 2-flux models. Figure 4(b) shows the equivalent q of the 2-flux model calculated using IBA a values. A comparison between this equivalent q and a demonstrates the complex relationship between them. At the snow surface, the equivalent q is very close to the value of a for all three pec, which according to the IBA, is 0.501, 0.514 and 0.543. However, it gradually increases from the snow surface to the bottom of 30 the snow profile. The equivalent q can reach larger values if the snow thickness is larger or the snow particles are larger. It is consistent with the conclusions in Figure 3 that, a strongly-forward scattering phase function is a sufficient condition to ensure low errors in approximating the 2-flux theory by the 1-flux theory. Of course, if q, the scattered ratio of intensities, is only considered as a parameter to describe the volume scattering event, it is reasonable to seek a value of q that gives the same simulated TB as the 2-flux model for a certain extent of snow micro-structure. However, it is true that, q is not as ideal as a, the forward-scattering ratio of phase function, because it is not constant along the snow profile, and is influenced by many factors, including snow depth, κs and ground and sky boundary conditions: see equation (2.11). Figure 5 shows the difference in TB at the snow surface between the 2-flux model and the 1-flux model with the default q (as 0.96) at 37 GHz. In Figure 5(a), results were shown versus snow slant thickness, d’. When snow slant thickness is 20 cm, the TB difference between two models at the snow surface is 18.3 K; however, when snow slant thickness is 2 meters, the difference is 71.9 K. After the first 2 meters, the decrease of TB in the 2flux model begins to saturate, but the TB of the 1-flux model continues to decrease. Figure 5(b) shows that the model difference increases with the increase of pec. In summary, generally in all cases, the simulated TB of 2-flux model is smaller than the 1flux model, because the forward scattered ratio of intensities, q, is fixed to 0.96, like in the HUT model, whereas the value of a is much closer to 0.5, as listed in Table 1. 31 250 200 150 100 TBup at the snow surface (K) TBup at the snow surface (K) 300 50 0 (a) 1−flux 2−flux 100 200 300 400 500 Snow slant thickness, d’ (cm) 600 700 800 300 250 200 150 100 80 0.05 (b) 1−flux 2−flux 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 Exponential correlation length, p (mm) ec 0.3 0.32 Figure 5. Upward TB at the snow surface simulated by the 1-flux and 2-flux models using the same absorption and scattering coefficients by Improved Born Approximation at 37 GHz: (a) TB versus snow slant thickness for medium snow grain size (pec=0.18 mm), (b) TB versus exponential correlation length (pec) for 50 cm snow slant thickness. 2) Scattering Coefficient Figure 6 shows the scattering coefficient for pec from 0.05 to 0.32 mm and Dmax from 0.24 mm to 4.74 mm at 37 GHz used by the HUT model and MEMLS. In general, the two HUT methods produce much larger scattering coefficients than the two MEMLS methods. The κs from the Hallikainen equation of the HUT model is 5.4 to 86 times of the κs from IBA, and is 2.6 to 7.8 times of the κs from empirical MEMLS. The κs from the Roy equation is 3.14 to 139.5 times of the κs from the IBA, and is 1.5 to 12.5 times of the κs from the empirical MEMLS. The extreme large ratios were found by comparing IBA with the HUT methods for pec at 0.05 mm, where IBA produced quite small κs with 32 reasons remaining to be studied. Such large contrast between the HUT and MEMLS κs was not found in Tedesco & Kim (2006a), because in that paper, pec was calculated as half of Dmax, which is valid if the snow particles are closely-packed small spheres. According to equation (2.33), the ratio of pec to Dmax ranges from 0.1 to 0.25, which is significantly smaller than 0.5. 0.23 25 κs (1/m) 20 15 0.31 Snow grain size, Dmax (mm) 0.41 0.54 0.72 0.94 1.25 1.65 2.2 2.8 3.8 4.7 Hallikainen eq. Roy eq. IBA of MEMLS empirical MEMLS 10 5 0 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 Exponential correlation length, pec (mm) 0.3 0.32 Figure 6. Scattering coefficient used by the HUT model and MEMLS as a function of pec at 37 GHz. To explore the reason for large difference between HUT κs and MEMLS κs, studies on previous papers are needed where their methods were first built to calculate or measure this parameter. The HUT model was derived using active measurement of the extinction coefficient of snow slabs by a free-space system (Hallikainen et al., 1987). MEMLS was based on passive microwave experiments, where the backward-scattering coefficient was solved by TB observed from snow slabs over different underlying conditions, and the scattering coefficient was calculated assuming symmetric forward and backward 33 scattering (Wiesmann et al., 1998). A possible reason for the different κs between the HUT model and empirical MEMLS is due to the experiment setup: the experiment based on the passive approach is insensitive to scattering in forward direction, and it could largely underestimate the retrieved κs if it assumes a=0.5. The difference between the HUT κs and the κs from MEMLS-IBA may also be that the experiment based on the active approach is not only affected by scattering due to snow grains, but also by the diffraction arising from large structures inside the slabs which could cause a strong forward scattering peak. However, the fact that the 2-flux κsforward is integrated over the forward half space can significantly average out the scattering peak near 0° scattering angle. This is more likely related to an angular resolution issue. According to deltaEddington approximation (Joseph et al., 1976) and the simulation experiment using Mie theory for a large sphere (size parameter as 4000) (Matzler, 2004a), the forward scattering peak could be truncated without influencing the RT equation resolution in simulating radiation fluxes. In other words, a large κs with a large forward-scattering ratio could be equivalent to a smaller κs and smaller forward-scattering ratio through reasonable transformation. In this case, the use of the large κs is thus consistent with the use of the large q. The difference between the scattering in the HUT model and the scattering in the MEMLS may be not that significant. Figure 7 compares the backward-scattering coefficient of MEMLS with 4% (1-0.96) of the HUT scattering coefficient. The later could be considered as the equivalent backwardscattering coefficient of the HUT model. Simulations show that this parameter of the two models is quite closer compared to the large contrasts in scattering coefficient. When pec 34 is larger than 0.225 mm, the IBA backward-scattering coefficient even becomes larger than the two kinds of HUT backward-scattering coefficients. Therefore, the parameterization of κe and the parameterization of q could be highly-related in the HUT model. In other words, the highly forward-scattered q could be compensating the very Scattering coefficeint in backward direction (1/m) large κs of the HUT model in some sense. 0.23 5 4 3 0.31 Snow grain size, Dmax (mm) 0.41 0.54 0.72 0.94 1.25 1.65 2.2 2.8 3.8 4.7 Hallikainen eq. Roy eq. IBA of MEMLS empirical MEMLS 2 1 0 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 Exponential correlation length, pec (mm) 0.3 0.32 Figure 7. Backward-scattering coefficient of MEMLS and 4% of the scattering coefficient used by the HUT model as a function of pec at 37 GHz. Figure 8 shows snow surface TB simulated by the 1-flux using the HUT κs parameterizations (blue lines), and by the 2-flux model using κs and a of MEMLS calculated by the IBA (red lines) at 37 GHz. Here the 6-flux theory is not yet included in blue lines of the 2-flux model simulations. Compared to Figure 5(a), the large HUT κs strongly decreases the TB of 1-flux model simulation; and thus, reduces the difference between the 1-flux model using q as 0.96 and the 2-flux model. For example, at 40 cm slant thickness when using the Hallikainen equation for HUT and IBA for MEMLS, the 35 TB difference decreases to -13.9 K from 32.2 K; at 1 meter slant thickness of snow, it decreases to -35.1 K from 57.2 K (these values were also summarized in Table 3). However, at 2-meter slant thickness, the value of difference is changed to -60.8 K from 71.9 K, thus the absolute difference is not reduced a lot. Comparing the average of the two methods of the HUT model and the average of the two models of MEMLS, it shows they are close for snow slant thickness smaller than 2 meters, but the difference increases the snow becomes deeper. Thus, when the large κs from the HUT model is used together with the large q, the simulation of 1-flux and 2-flux models could still not closely-match. An interesting thing is that, at pec=0.18 mm when IBA and two methods of the HUT model give almost the same backward-scattering coefficient (Figure 7), the simulated TB at snow surface by IBA is still higher for slant snow depth over 80 cm (Figure 8(a)). The reason is that, inside a homogenous snow medium, the upward flux in the 1-flux model is not linked to downward flux; however, in the 2-flux model, the upward flux, as it propagates through snowpack, has contribution from reflected downward flux at all depth. This effect is stronger for deep snow, because in this case the downward flux could increase to higher values at snow bottom. Such mechanism stops the upward flux of the 2-flux model from dropping too low for deep snow. This is also why using close values of backward scattering coefficients, the saturation of the 2-flux model is still faster than the 1-flux model. 36 TBup at the snow surface (K) TBup at the snow surface (K) 300 (a) 250 200 150 100 50 0 100 200 300 400 500 Snow slant thickness, d’ (cm) 600 700 300 (b) 250 200 150 800 Hallikainen eq. Roy eq. IBA (without 6−flux) emp. MEMLS (without 6−flux) IBA (with 6−flux) emp. MEMLS (with 6−flux) 100 80 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 Exponential correlation length, pec (mm) 0.3 0.32 Figure 8. Upward TB at the snow surface at 37 GHz simulated by the 1-flux and 2-flux models using the different scattering coefficients of the HUT model (red lines) and MEMLS, with 6-flux theory included (green lines) and not included (blue lines): (a) TB versus snow slant thickness for medium snow grain size (pec=0.18 mm, correspondingly Dmax =1 mm), (b) TB versus exponential correlation length for 50 cm snow thickness. 3) Trapped Radiation The green lines in Figure 8 show the 2-flux simulations with 6-flux theory employed. It shows the consideration of trapped radiation in MEMLS increases its simulated TB. Using the IBA, the increase of TB due to the 6-flux theory is 10.4 K at 40 cm slant thickness, 21.6 K at 1-meter slant thickness, and 32.4 K at 2-meter slant thickness. The influence of the trapped radiation increases with the increasing snow depth, until it reaches a large depth of about 2-meter slant thickness. As in Table 2, for cases where the combined results of the first two factors (the different κs and the 1-flux versus 2-flux 37 theories) already gives higher MEMLS TB than HUT TB, considering trapped radiation in MEMLS increases the absolute TB difference between the models. For example, in average of four pairs of model approach comparisons, the TB difference was 86.8 K for 2 meter snow slant thickness if the same κs as that of MEMLS was used; after employing different κs, it was reduced to -19.5 K, with the 2-flux model slightly higher than the 1flux model; however, after including 6-flux theory, the TB difference changed to -57.2 K, with its absolute value increased. It means the predicted effect of trapped radiation is not negligible for deep snow. 4) Combined Effects Table 2 lists more values of TB differences between the HUT and MEMLS approaches at 37 GHz. In addition to the effects already mentioned above, it shows when all three factors are combined, in most cases, the simulated TB at the snow surface is higher using MEMLS than using the HUT model. The intra-differences between different κs parameterizations of the same model are also quite large. In the summarized averages (see the last 3 lines), large difference between complete HUT and MEMLS models is found (-57.2 K) for 2-m slant thickness snow, i.e. 1.15 m snow depth given 55º propagation angle. As shown in the red (HUT) and green (6-flux MEMLS) lines in Figure 8, it is due to the summed effects of all three model theoretical differences mentioned above. After understanding how three factors work to give the final HUT and MEMLS TB results, it is time to address the difference in performance of full-version models, as will be shown in next section. 38 Difference of Simulated TB above the Snow Surface between models (K) TB(Halli.)– 40-cm thickness TB(IBA) 1-m thickness 2-m thickness TB(Roy) – 40-cm thickness TB(IBA) 1-m thickness 2-m thickness TB(Halli.)– 40-cm thickness TB(emp.) 1-m thickness 2-m thickness TB(Roy) – 40-cm thickness TB(emp.) 1-m thickness 2-m thickness Average 40-cm thickness 1-m thickness 2-m thickness Using the same κs 32.2 57.2 71.9 32.2 57.2 71.9 70.9 100.1 101.7 70.9 100.1 101.7 51.6 78.7 86.8 Using different κs of HUT and MEMLS parameterizations -13.9 -35.1 -60.8 -0.8 -11.5 -32.2 31.2 22.1 -6.8 44.3 45.7 21.8 15.2 5.3 -19.5 Using different κs and including 6-flux theory -24.3 -56.7 -93.2 -11.2 -33.1 -64.6 15.5 -9.6 -49.8 28.6 14 -21.2 2.2 -21.4 -57.2 Table 2. Difference of simulated TB above the snow surface for 40-cm, 1-meter and 2meter slant snow thickness with medium snow grain size at 37 GHz. Section 2.3: Comparison of Model Simulations with Passive TB Measurements The numerical experiments in previous sections focused on isolating the implications of three model differences for idealized snow, but truly the snow cover in the field is never ideal. Therefore, here the performance of the HUT model and MEMLS was compared against in-situ observations. The snow datasets used for comparison include: (1) measurement at the Intensive Observation Area (IOA) of the first and second Nordic Snow Radar Experiment (NoSREx) carried out at Sodankyla, Finland (67.362°N, 26.633°E) during the winters 2009, 2010, and the follow-on experiments of 2011 and 2012 (Lemmetyinen et al., 2016); (2) measurement of 8 snowpits in the clearings of a forest region near Churchill, Manitoba, Canada (58.73°N, -93.82°E) from November 2009 to April 2010 (see forest sites in Derksen et al. (2012)); (3) measurement of 6 snowpits during the Intensive 39 Observation Period 3 (IOP3) of the Cold Land Processes Field Experiment (CPLX) at the Local Scale Observation Site (LSOS) (39.9066°N, -105.8829°E) near Fraser, Colorado in February 19-24, 2003 (Hardy et al., 2003); (4) measurement of a 174-cm snowpit at Storm Peak Laboratory (SPL), Colorado (40.455°N, -106.744°E) in February, 2010 (Durand et al., 2012). In the snow classification system (Sturm et al., 1995), the snow at Colorado is alpine snow, where the snow at Sodankylä and Churchill are taiga snow. In these experiments, the snow layer thickness, temperature and density were all measured. For the snow microstructure, all snowpits have Dmax measurements. Besides Dmax, the snowpit at SPL, Colorado has pec measured by exponentially fitting the spatial autocorrelation function using the stereology image (Matzler, 2002). The dataset at Sodankyla has two snowpits with SSA measured by micro-CT (Computer Tomography) (Schneebeli et al., 2004). The measured SSA was converted to pec using the methods in Matzler (1997) and (2002). For the condition of the underlying soil, the soil physical temperature or the temperature at the snow-soil interface was measured at all sites. The soil moisture was measured via the TDR sensor at 2 cm below soil surface at Sodankyla, and at 1.5 cm below soil surface at LSOS in CLPX. Above the snow surface, ground-based radiometers were used to measure the TB from snowpack. Table 3 lists the details pertaining to the radiometric observations. 40 Site Radiometer Frequency (GHz) Observation angle Mounted height 3-dB beam width Sodankyla (Lemmetyinen et al., 2016) Sodrad (FMI Sodankylä Radiometer) 10.65, 18.7, 36.5, 90 50° 5 m above the soil surface < 6.1° 6.9, 19, 37, 89 53° 1.55 m above the snow surface 6-9° 18.7, 89 55° 5.8 m above the soil surface 10° 50° 3.5 m above the soil surface 8° Churchill (Derksen et al., 2012) CLPX (Graf et al., 2003) SPL Environment Canada sledbased microwave radiometer system GBMR (Ground Based Microwave Radiometer)-7 AESMIR (Airborne Earth Science Imaging Radiometer) 36.5, 18, 36 Table 3. Details of TB measurements at all sites. For the simulations of the HUT model, all the measured snow grain size (Dmax) was processed into the effective snow grain size (Deff) using equation (2.29). For the simulation of MEMLS, the pec measured at SPL by stereology was used directly; at other sites, equation (2.33) was used to convert Dmax into pec. In equation (2.33), the coefficients fitted from Weissflühjoch snow were used for CLPX snowpits at Colorado, and the Churchill snowpits. However, at Sodankylä, it showed that, the observed pec were slightly larger than the predicted values using the Weissflühjoch coefficients. Especially, because the snow density of taiga snow is typically lower than alpine snow (Sturm et al., 1995), the density threshold v0, smaller than which the snow has a relatively stable p0 that is unrelated to Dmax, is smaller at Sodankylä than at Weissflühjoch. Using the original coefficients result in overestimation of TB on several days when the snowpits indicated a 41 low-density surface layer. Therefore, a new set of coefficients were refitted for use at Sodankyla, as: a0 = 0.23 mm, a1 = 0.13 mm, p0 = 0.09 mm, v0 = 0.16 and d0 = 0.5 mm. The fitted errors are: εfit = 0.032 mm, and ε0 = 0.012 mm. “Ice layers” were recorded in the field for some CLPX snowpits, and some snowpits at Churchill, but described differently at different sites. In this paper, previous studies using the same datasets were referred to for the treatment of ice layers. As in Durand et al. (2008), the “ice layers” at CLPX were considered as melt-refreeze crusts, with pec as 0.25 mm, Dmax as 2.2 mm, and snow density as 500 kg m-3. The “ice layers” at Churchill were considered as pure ice with pec and Dmax as 0 mm, and density as 917 kg m-3 as in Montpetit et al. (2013). The air-snow boundary and the interfaces between snow layers were all assumed smooth, and the reflectivities at these surfaces were simulated using Fresnel equations, following common practice (Lemmetyinen et al., 2010). The reflectivity of soil-snow boundary was calculated referring to previous studies using the same datasets. A standard atmosphere was used to simulate to obtain the downward sky TB at each frequency and incidence angle as described by Durand & Margulis (2007). For example, at 55° incidence angle, Tsky is 5 K, 6 K, 21 K, 33 K and 94 K at 6.925, 10.65, 18.7, 36.5 and 89 GHz, respectively. The snowpits in forest clearings may be influenced by the surrounding forest canopy, which emits downward radiation upon the surface of snowpits in addition to the sky emission. If the vegetation plus sky emission reach 51 K at 19 GHz, and 62 K at 37 GHz according to Montpetit et al. (2013), the influence to the HUT simulation is less than 2 K; for MEMLS, the effect is higher (3-11 K at 37 GHz, 3-7 K at 42 19 GHz) due to the stronger relation between the upward and downward fluxes in the 2flux theory. These values are given for reference, as the focus of the study is to understand the different treatment of microwave scattering within the snow rather than model validation. Without detailed information regarding tree height, temperature, stem volume, distance to the snowpits, etc., the effect of vegetation cannot be thoroughly addressed, and thus is not considered further in this paper. After all simulations were completed, the errors between modeled and observed TB were analyzed. Figure 9 shows the simulated versus observed TB at 37 GHz, vertical polarization across all sites and all years. Four sets of simulation approaches were used: the HUT model using κs from the Hallikainen (a) and Roy (b) equations, and MEMLS using both IBA (c) and empirical (d) κs. First, note that when both MEMLS and the HUT model show some skill in reproducing the observed TB, the best of the model cases has an R2 of 0.46, explaining less than half of the variance in observed microwave brightness. One major reason for the lack of skill compared with other studies (e.g. Liang et al. (2008), Rutter et al. (2014)) is that we have not applied any calibration to the grain size to better fit the TB observations. Also, most of the simulations shown here used geometric grain size, which has uncertain links to the effective scatterer sizes. Statistics, however, show that the RMS error for the MEMLS using IBA case is 11.0 K at 37 GHz, vertical polarization, which is close to that found in past studies (Durand et al., 2008; Montpetit et al., 2013). Second, note that MEMLS generally has higher R2 values than HUT, despite the possible introduction of uncertainty by converting Dmax to pec. Thirdly, it shows that the MEMLS using IBA has a lower TB bias than MEMLS using the empirical scattering 43 coefficient. Similarly, the HUT model using the Roy equation has a lower TB bias than the HUT model using the Hallikainen equation. The models all underestimated TB at Sodankylä for the third winter season. The field data indicated thicker snow depth and smaller bottom snow grain size in that winter, but these features still could not reproduce the even higher observed TB. At CLPX, the experiment was taken basically within one week. Although there was new snowfall in the middle of that week (Durand et al., 2008), the scattering of surface new snow layer is not strong. Thus, the variation in observed TB at 37 GHz, vertical pol. was within 3.5 K, which makes it difficult for the models to capture its features. The rms error using MEMLS-IBA at CLPX is 5.6 K at 37 GHz. 44 (a) HUT(Halli.) 260 R2 =0.289 240 reg Simulated TB(K) Simulated TB(K) 240 260 220 200 180 160 260 Simulated TB(K) 240 160 180 200 220 Observed TB(K) 240 220 200 180 140 140 260 (c) MEMLS(IBA) 260 R2 =0.525 reg 2 240 R =0.459 220 200 180 160 140 140 R2reg=0.237 160 Simulated TB(K) 140 140 (b) HUT(Roy) 220 200 160 180 200 220 Observed TB(K) 240 260 240 260 (d) MEMLS(emp.) Sdkl1 Sdkl2 Sdkl3 Sdkl4 Churc CLPX SPL R2 =0.438 reg 180 160 160 180 200 220 Observed T (K) B 240 260 140 140 160 180 200 220 Observed T (K) B Figure 9. The simulated TB of MEMLS and HUT compared to TB observations, at 37 GHz in vertical polarization, where the black dash line is the 1:1 line. The R2 is the coefficient of determinant between the simulation and the observation (i.e., the R2 of the 1:1 line), and the R2reg is the coefficient of determinant for the fitted lines between the observed and the simulated TB. To study the difference of the HUT model and MEMLS for deep snow, Figure 10 plots the simulated bias in TB as a function of snow depth. As in Fig. 10(a) and (b), snow depth empirically explains a significant amount of the variability in TB simulation errors for the 45 HUT model. It shows that, the HUT model (both κe parameterization) produces an increasingly underestimation of TB at 37 GHz as snow depth increases; where the MEMLS errors (both ways of scattering coefficient calculation also) are unrelated to snow depth. If using linear relationship to fit TB bias, dTB, as a function of snow depth, SD, purely for the purpose of data analysis, it could be written as: dTB [K] = a * SD [cm]+ b (2.34) At 37 GHz, vertical polarization, the coefficients for Hallikainen equation are: a = -0.662 K/cm, b = 14.14 K, with R2=0.551; the coefficients for Roy equations are: a = -0.548 K/cm, b = 28.42 K, with R2=0.550. According to these equations, if acceptable absolute TB error is 20 K, the underestimation needs to be considered for snow depth larger than 52 cm for the Hallikainen equation, and 88 cm for the Roy equation; if it is 10 K, the snow depth needs to be no larger than 36 cm for the Hallikainen equation, and 70 cm for the Roy equation. 46 −20 −20 −40 −60 −60 −80 −80 20 B dTB(K) 0 −40 50 100 150 Snow depth (cm) −100 0 200 (c) MEMLS(IBA) 20 0 0 −20 −20 −40 −60 −80 −80 50 100 150 Snow depth (cm) 200 50 100 150 Snow depth (cm) −100 0 200 (d) MEMLS(emp.) −40 −60 −100 0 (b) HUT(Roy) 20 0 −100 0 dT (K) (a) HUT(Halli.) dTB(K) B dT (K) 20 Sdkl1 Sdkl2 Sdkl3 Sdkl4 Churc CLPX SPL 50 100 150 Snow depth (cm) 200 Figure 10. Bias of MEMLS and HUT simulations compared to TB observed by radiometers as a function of snow depth at 37 GHz in vertical polarization. The black dash line is where dTB equals zero. Figure 11 shows that the underestimation of HUT at 37 GHz in Figure 10 is uncorrelated to snow grain size. Of all the methods, only empirical MEMLS shows some decreasing trend of dTB to snow grain size, but with insignificant R2 of 0.21. 47 (a) HUT(Halli.) 0 0 −20 −20 −40 −40 −60 −60 −80 −80 −100 0 20 0.5 1 1.5 Snow grain size (mm) −100 0 2 (c) MEMLS(IBA) 20 0 −40 1 1.5 2 (d) MEMLS(emp.) −20 Sdkl1 Sdkl2 Sdkl3 −60 Sdkl4 −80 CLPX −100 0 0.5 Snow grain size (mm) 0 dTB(K) B dT (K) −20 (b) HUT(Roy) 20 dTB(K) B dT (K) 20 −40 −60 Churc −80 SPL 0.5 1 1.5 Snow grain size (mm) 2 −100 0 0.5 1 1.5 Snow grain size (mm) 2 Figure 11. Bias of MEMLS and HUT simulations compared to TB observed by radiometers as a function of snow grain size at 37 GHz in vertical polarization. The black dash line is where dTB equals zero. Figure 12 shows the sensitivity of TB at 19 GHz in vertical polarization to snow depth. At this band, the HUT model using the Hallikainen equation shows almost no underestimation for deep snow, whereas the significant decreasing trend remains the 48 same as at 37 GHz with the Roy equation. Combining Figure 10 and Figure 12, the TB difference between 19 and 37 GHz increases with snow depth if the Hallikainen equation is applied to both frequencies, because 37 GHz is biased low for deep snow whereas 19 GHz is not. Using the Roy equation, a bias of about -30 K appears for all depth in the TB difference between 19 and 37 GHz, but a good thing is that its dependency to snow depth cancels out. For sensitivity analysis of snow grain size (Figure 13), the conclusions at 19 GHz are close to those at 37 GHz: the HUT model and the MEMLS-IBA performances are insensitive to snow grain size. 49 (a) HUT(Halli.) 0 0 −20 −20 −40 −40 −60 −60 −80 −80 −100 0 20 50 100 150 Snow depth (cm) −100 0 200 (c) MEMLS(IBA) 20 0 Sdkl1 100 150 Snow depth (cm) 200 (d) MEMLS(emp.) Sdkl2 −40 Sdkl3 −60 Sdkl4 −80 CLPX −20 dTB(K) B dT (K) 50 0 −20 −100 0 (b) HUT(Roy) 20 dTB(K) B dT (K) 20 −40 −60 Churc −80 SPL 50 100 150 Snow depth (cm) 200 −100 0 50 100 150 Snow depth (cm) 200 Figure 12. Bias of MEMLS and HUT simulations compared to TB observed by radiometers as a function of snow depth at 19 GHz in vertical polarization. The black dash line is where dTB equals zero. 50 (a) HUT(Halli.) 0 0 −20 −20 −40 −40 −60 −60 −80 −80 −100 0 20 0.5 1 1.5 Snow grain size (mm) −100 0 2 (c) MEMLS(IBA) 20 0 −40 1 1.5 2 (d) MEMLS(emp.) −20 Sdkl1 Sdkl2 Sdkl3 −60 Sdkl4 −80 CLPX −100 0 0.5 Snow grain size (mm) 0 dTB(K) B dT (K) −20 (b) HUT(Roy) 20 dTB(K) B dT (K) 20 −40 −60 Churc −80 SPL 0.5 1 1.5 Snow grain size (mm) 2 −100 0 0.5 1 1.5 Snow grain size (mm) 2 Figure 13. Bias of MEMLS and HUT simulations compared to TB observed by radiometers as a function of snow grain size at 19 GHz in vertical polarization. The black dash line is where dTB equals zero. Table 4 and Table 5 summarize the mean bias and rms error, across all sites, and for the various models and scattering coefficient parameterizations, and across both polarizations and all frequencies. They show that, for the two methods of HUT, the Hallikainen equation is better for low frequencies, including 6.925, 10.65 and 18.7 GHz; the Roy 51 equation is better for high frequencies, including 36.5 and 90 GHz. Especially at 90 GHz, the underestimation of the Hallikainen equation reached 100 K. In Hallikainen et al. (1987), this κe parameterization at 90 GHz showed large divergence compared to the theoretical values by strong fluctuation theory, too. Of the two methods for MEMLS, the IBA appears to have lower errors, where the empirical MEMLS usually shows underestimation. The performance of the models is better at vertical polarization than horizontal polarization, because horizontal polarization is often impacted by boundary reflectivities and sub-layer inhomogeneity of snow properties, etc. Improved treatment of ice layers is important for model improvement. For example, considering ice layers according to Montpetit et al. (2013) at Churchill reduces the error in estimating polarization difference at 19 GHz; at Sodankylä, underestimation of polarization difference was found at all frequencies, but it could be partly reproduced where ice layers were reported. Bands V7 V11 V19 V37 V89 H7 H11 H19 H37 H89 HUT(Halli.) 6.0 5.6 -4.6 -26.9 -111.7 -29.0 10.5 5.5 -16.8 -101.7 Mean Bias (K) MEMLS HUT(Roy) (IBA) -14.8 6.0 -18.9 7.7 -31.6 1.8 -5.6 -1.4 -4.6 2.8 -44.4 -35.2 -13.2 12.7 -19.6 12.4 3.4 9.7 1.5 11.3 MEMLS (emp.) 1.8 -6.9 -33.8 -31.4 -3.5 -37.1 -0.5 -19.2 -18.3 5.3 Table 4. Mean biases of the simulated compared to the observed TB above the snow surface at all sites 52 Rms Error (K) Bands V7 V11 V19 V37 V89 H7 H11 H19 H37 H89 HUT(Halli.) HUT(Roy) 7.1 6.6 9.2 33.9 117.3 38.4 13.8 17.2 27.7 106.9 16.8 20.9 35.0 17.9 29.5 53.7 17.7 27.6 19.2 27.0 MEMLS (IBA) 7.6 8.4 6.9 11.0 17.1 47.1 15.5 19.1 16 19.9 MEMLS (emp.) 5.1 10.1 37.4 34.4 19.3 48.9 10.6 26.4 23.0 18.8 Table 5. Rms errors of simulated compared to the observed TB above the snow surface at all sites In summary, results of this section indicate that HUT can be expected to predict much lower TB than MEMLS for deep snow. The reason for this is that even if the scattering coefficient and phase function between the HUT model and MEMLS can be shown to be equivalent, the error in the 1-flux simplification and the trapped radiation are still nonnegligible for thick snow. Comparison with in-situ observations at 37 GHz (vertical polarization) shows that MEMLS using the κs from Improved Born approximation (IBA) has the best performance, with RMSE as 11.0 K, and R2 as 0.46. No calibration or tuning of the grain size was performed. MEMLS using empirical parameterization of κs underestimates the TB at the snow surface. The HUT model tends to underestimate TB at the snow surface for thick snow, and there is strong correlation between the TB bias and snow depth, with R2 over 0.5 for its both κs parameterizations. Therefore, in the 4 option of snow RT models, MEMLS-IBA is more accurate than other models, and is suggested to be used for SWE retrieval based on passive microwave measurements. 53 Section 2.4: Comparison of Model Simulations with Active Backscattering Measurements The MEMLS model has been adapted to include calculation of the backscattering coefficient calculation, as published in Proksch et al. (2015). In the previous section, the Improved Born Approximation (IBA) performed better than the empirical MEMLS for volume scattering coefficient calculation. Therefore, here MEMLS3&a is used with IBA to simulate backscattering coefficients for the Sodankyla snowpits (Lemmentynien et a., 2016). During the Nordic Snow Radar Expeirments (NoSREx), the active backscattering coefficients at three frequencies (10.2, 13.3 & 16.7 GHz), four polarizations (VV, HH, VH & HV) were observed by the SnowScat (Snow Scattermoeter) for the same snowpit measurements utilized in the previous section. The observation angles include 30 to 60˚ per 10˚ step. The model comparison here focused on 50˚. MEMLS3&a Introduction The main assumption of MEMLS3&a is the bistatic scattering coefficient at directions other than the forward-scattering direction is diffuse and Lambertian (Proksch et al., 2015): σ d0 = 4rd u02 (2.35) where u0 = cosθ 0 is the cosine of the radar incidence angle (also the observation angle in the backscattering direction), rd represents the diffuse component of the total reflectivity observed at the snow surface, r, and rd = r − rs , where rs is the specular component of the total snowpack reflectivity. Accordingly, σ d is the backscattering coefficient component 54 resulted from diffuse scattering in the snowpack, which may include the diffuse snow volume scattering, the diffuse part of air-snow, snow-soil surface scattering and their interactions with the snow volume scattering. Note that σ d includes both the copolarized (VV, HH pol.) and cross-polarized (VH, HV pol.) backscattering signals. MEMLS3&a uses an empirical splitting parameter q to calculate σ d,0 pp' at pp’ polarization as: ⎧ 0 (1− q)σ d,v ⎪ ⎪ 0 σ d,0 pp' = ⎨ (1− q)σ d.h ⎪ ⎪⎩ q(σ d,v + σ d,h ) / 2 p=p'=v p=p'=v (2.36) p=v, p'=h; or p=h, p'=v 0 0 where, σ d,v and σ d,h are calculated by rdv and rdh using equation (2.35) for each specific polarization. There is a close relationship between MEMLS and MEMLS3&a, such as the total reflectivity at snow surface, r, is calculated from the TB simulated by MEMLS at the same frequency as (Proksch et al., 2015): r= TB1 − TB2 Tsky1 − Tsky2 (2.37) where, TB1 and TB1 are the brightness temperature observed at the snow surface simulated by exactly the passive MEMLS, but they used different downwelling sky brightness temperature (i.e. Tdn(0+)) in order to calculate r. TB1 uses Tsky1 (=100 K) and TB2 uses Tsky2 (=0 K). The next step is to calculate the specular part of snow surface reflectivity rs. As in Proksch et al. (2015), MEMLS3&a calculates rs as the specular part of soil surface 55 reflectivity attenuated by snow layers. Snow layer interfaces were assumed smooth, and multiple reflections between layer interfaces were considered by a recurrence approach. See the Appendix in Proksch et al. (2015) for details. The attenuation factor of each snow layer is calculated as the sum of the MEMLS absorption coefficient and 4 times of the MEMLS scattering coefficient; the later contains an empirical scalar of 4 to consider the diffraction effect of coherent waves. Some part of the specular reflectivity at the snow surface can be scattered to the backscattering direction when the snow surface has surface undulations (roughness). The specular part of backscattering coefficient, σ s0 , is calculated using Geometrical-Optics (GO) theory for small roughness in MEMLS3&a as: 0 s σ = rs,0 exp (−tan 2 θ 0 / (2m 2 )) 2m 2 µ 04 (2.38) where m2 is the mean-squared slope (rms-height/correlation length) of surface roughness, θ 0 and µ0 are related to the backscattering angle, rs,0 is the specular reflectivity at normal incidence calculated as the mean of rs,0v and rs,0h at vertical and horizontal polarizations. In the end, the total backscattering coefficient measured by the radar is simulated as: σ 0 pp' ⎧ 0 0 ⎪ σ d, pp' + σ s =⎨ 0 σ ⎪ ⎩ d, pp' p=p' (2.39) p ≠ p' In summary, several simplifications were made in MEMLS3&a: (1) The incident radar pulses have narrow beamwidth, which physically requires fineresolution RT modeling with a large number of fluxes. However, by assuming Lambertian diffuse scattering, MEMLS3&a adapted the passive MEMLS, which 56 essentially is a 6-flux model, to simplify the RT equation resolution. The assumption of Lambertian diffuse scattering may be violated by the backscattering enhancement mentioned in Tan et al. (2015), which happens for rough soil surfaces due to the coherency of double-bounce waves in the backscattering direction. (2) MEMLS3&a uses an empirical factor q to calculate the cross-pol. backscattering coefficient. In reality, this effect is physically be related to the non-spherical shape of snow particles. Note that the q here is not the forward-scattered ratio of intensities used in HUT model. Hereafter, a new notation qs is used for the polarization splitting factor to prevent ambiguity. (3) MEMLS3&a uses m to calculate the backscattered specular radiation. However, because r and rs are calculated for the entire snowpit, m more likely also represents the overall undulation of the entire snowpit, which is not limited to the snow surface roughness, but also includes the roughness of snow layer interfaces and sub-layer structures. MEMLS3&a Simulation Results Figure 14 shows the performance of MEMLS3&a to simulate the backscattering coefficient at VV-pol. at Sodankyla for 4 Intensive Observation Periods (IOP) from 2009 to 2013. The simulations here used the measured multi-layer snow properties and the measured soil temperature and moisture. The soil surface roughness was fixed as 0.1 cm, because the soil surface was reported smooth (Lemmetyinen et al. 2016). The MEMLS3&a model parameters are: m as 0.1 and qs as 0.1 for all frequencies. Simulation 57 tests show that m only influences backscattering coefficient measured at less than 40˚ observation angle. If pec was scaled to match the measured TB at 36.5 GHz, vertical polarization for each snowpit, a qs between 0.8 and 1.2 can be calculated for different snowpits and frequencies. However, such simulations are based on a snowpit-specific optimization for pec, and do not consider other factors which may also influence backscattering coefficient (eg. soil conditions, snow sub-layer structures). Therefore, here an intermediate value of qs = 0.1 is used. First, it should be noted that no specific scaling of pec were applied for different snowpits in Figure 14. Therefore, not all the simulations at three frequencies match the observations. Secondly, in the early season of all IOPs except IOP2, the simulated backscattering coefficient at three frequencies is smaller than the observation, with reasons to be explained later with Figure 16 to Figure 19. There are two assumptions to explain these underestimations: one is related to the possible large snow grain size at the bottom of the snowpit, developed from the melt-refreezing events at the early snow season; the other is related to the soil background emission, because the underestimation is found also at 10.2 GHz. Thirdly, it can be found that MEMLS3&a underestimates 16.7 GHz VV-pol. backscattering coefficient for the entire IOP2, which is not found for other IOPs. Such underestimation can be fixed using a larger pec. Therefore, the reason for the underestimation at 16.7 GHz is more likely that scattering from large crusts occurred in the snowpacks during IOP2, or that there were uncertainties in converting the measured Dmax to the MEMLS pec. IOP2 has the largest snow grain size (Dmax) observed at the bottom of the snowpacks among all IOPs. For large snow particles, a small error in pec 58 can result in a large error in snow volume scattering calculation. Fourthly, it should be noted that for IOP4, there is a discontinuous point for the observation at 10.2 GHz on February 28. This is because of a reinstallation of the SnowScat instrument at that time. IOP1 (2009-2010) - VV -8 -10 VV (dB) 0 -12 -14 -16 -18 -20 -22 Oct Nov Dec Jan Feb Mar Apr May Jun Jul Apr May Jun Jul IOP2 (2010-2011) - VV -8 -10 VV (dB) 0 -12 -14 -16 -18 -20 -22 Oct Nov Dec Jan Feb Mar Continued Figure 14. Observed SnowScat VV-pol. backscattering coefficients (lines) at Sodankla compared to MEMLS3&a simulations (crossed) at four Intensive Observation Periods (IOP). The blue, red and green colors are for 10.2, 13.3 and 16.7 GHz, respectively. 59 Figure 14 continued IOP3 (2011-2012) - VV -8 -10 VV (dB) 0 -12 -14 -16 -18 -20 -22 Oct Nov Dec Jan Feb Mar Apr May Jun IOP4 (2012-2013) - VV -8 10.2-Obsr 13.3-Obsr 16.7-Obsr 10.2-Simu 13.3-Simu 16.7-Simu -10 -12 VV (dB) 0 Jul -14 -16 -18 -20 -22 Oct Nov Dec Jan Feb Mar Apr May Jun Jul Figure 15 shows the MEMLS3&a simulations at VH polarization. First, it can be seen that most of the simulations match well with the observations. It indicates qs as 0.1 is suitable for the Sodankyla snowpits. Secondly, some underestimations are still found in early snow season for IOP1, 3 and 4, like VV-pol. simulation results in Figure 14. 60 IOP1 (2009-2010) - VH -15 VH 0 (dB) -20 -25 10.2-Obsr 13.3-Obsr 16.7-Obsr 10.2-Simu 13.3-Simu 16.7-Simu -30 -35 Oct Nov Dec Jan Feb Mar Apr May Jun Jul Apr May Jun Jul Apr May Jun Jul IOP2 (2010-2011) - VH -15 VH 0 (dB) -20 -25 -30 -35 Oct Nov Dec Jan Feb Mar IOP3 (2011-2012) - VH -15 VH 0 (dB) -20 -25 -30 -35 Oct Nov Dec Jan Feb Mar Continued Figure 15. Observed SnowScat VH-pol. backscattering coefficients at Sodankla compared to MEMLS3&a simuations at four IOPs. 61 Figure 15 continued IOP4 (2012-2013) - VH -15 VH 0 (dB) -20 -25 -30 -35 Oct Nov Dec Jan Feb Mar Apr May Jun Jul Table 6 summarizes the accuracy of MEMLS3&a simulations. Due to the underestimation at the early snow season mentioned above, simulated backscattering coefficient at VV-pol. is biased about -1.15 to 0.02 dB at 10.2 GHz, -0.91 to 0.4 dB at 13.3 GHz and -2.2 to -0.03 dB at 16.7 GHz for different IOPs. The rms error of VV-pol. backscattering coefficient is about 1 to 2.7 dB. The rms error of VH-pol. simulations is close to the accuracy of VV-pol. simulations. Similar to the TB results shown previously, the simulations are less accurate at HH-pol. due to the influence of horizontallyorientated sub-layer structures in snow. 62 Time IOP1 IOP2 IOP3 IOP4 Freq. (Pol.) Mean bias (dB) VV VH 10.2 -0.90 13.3 16.7 Rms error (dB) HH VV VH HH Number -0.03 2.48 1.43 1.19 2.50 -0.57 0.00 2.37 1.41 1.41 2.45 -1.67 -1.45 0.07 2.00 1.93 0.99 10.2 0.02 0.57 2.95 1.21 1.21 3.10 13.3 -0.33 0.41 2.64 1.09 1.13 2.74 16.7 -2.20 -1.01 -0.01 2.36 1.39 0.68 10.2 -1.15 -0.79 1.65 1.30 1.10 1.56 13.3 0.40 0.32 2.93 1.09 1.31 2.53 16.7 -0.03 0.29 2.18 0.95 1.20 2.04 10.2 -0.67 0.69 2.83 2.10 2.18 3.29 13.3 -0.91 0.23 2.13 2.70 2.69 3.08 16.7 -0.49 -0.12 1.75 2.73 2.89 3.08 30 20 10 23 Table 6. Mean bias and rms error of the MEMLS3&a simulations for the Sodankyla snowpits. Discussion To understand the MEMLS3&a simulation results, Figure 16 to Figure 20 show a series of in-situ snow and soil measurements and map them together with the SnowScat measurements. The measured frost soil depth (Rautiainen et al., 2014) and the soil moisture measured at 5- to 80-cm depths (Lemmetyinen et al., 2016) provide deep soil information where the radar pluses at 10.2 GHz may reach. The snowcover onset, snowmelt and the important changes of soil and air temperatures are indicated in these figures. Early in the snow season for all IOPs, the underlying soil gradually froze from top to bottom; at the same time, the backscattering coefficient decreased with slightly increasing snow depth. When the frost soil depth reached about 80-cm to 120-cm and the soil moisture at 2-cm to 80-cm became stable, the decreasing trend for backscattering 63 coefficient stopped, and turned to increase with increasing snow depth due to the increased amount of snow volume scattering. On the other hand, the measured snow grain size (Dmax) at the early snow season was no larger than the grain size in the middle of the snow season. Therefore, the higher measured SnowScat backscattering coefficient than the MEMLS3&a simulation at the early snow season is more likely related to the complex soil conditions beneath the snow. The partly frozen soil has a high-moisture layer beneath the surface frozen layer, which can produce strong surface scattering signals. Frozen soil may also have some extent volume scattering effect, which can increase the backscattering. In MEMLS3&a simulation, only the soil moisture at 2-cm was used to calculate backscattering coefficient. Based on an assumption of half-space homogenous soil medium, it didn’t essentially consider the layered effects inside soil. 64 IOP1 2009-2010 -20 Snowcover Onset 160 Tsoil2cm<0 o C 120 80 40 0 -40 -80 Feb Mar Apr May 100 50 0 -50 -120 -160 -200 Oct -100 Nov Dec Jan Feb Mar Apr May Jun Jul 1 0.1 Snow Melted 0.15 Mean Tair>0 o C Tsoil2cm<0 o C 0.2 SoilW-2cm FrostDepth>80cm Snowcover Onset 0.4 0.2 60 Mar Apr May 40 20 Nov Dec Jan Feb Mar Apr Jun Snow Melted 80 Feb Tair>0o C Tsoil2cm<0 o C 100 Snow Depth (cm) Jan Mean Dec Snowcover Onset Nov 120 0 Oct 0.8 0.6 0.05 FrostDepth>80cm Soil Water Content (m 3 /m 3 ) 0.3 0.25 0 Oct GHz (dB) 0 Jul Jun Snow Depth (cm) Jan 4 2 Snow Melted Dec Mean Tair>0 o C Nov FrostDepth>80cm -22 Oct 200 GHz - 10.2GHz 13.3GHz 16.7GHz 16.7GHz - 10.2GHz 6 VV 16.7 0 8 VV 10.2 0 10 -16 -18 Soil Frost Depth (cm) 12 Snow Melted -14 14 Mean Tair>0 o C VV (dB) 0 -12 FrostDepth>80cm -10 Tsoil2cm<0 o C Snowcover Onset -8 May (mm) Dmax 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 Jun 0 Jul 120 100 80 60 40 20 0 Jul Figure 16. Measured backscattering coefficient, snow depth & soil frost depth, soil moisture content and snow grain size (Dmax) (from top to bottom) at Sodankyla during IOP1. The gray vertical lines mark important events from snowcover, soil temperature and frost depth measurements; The red dash lines indicate a change of features for the backscattering coefficient measurement. 65 Snow Depth (cm) 0.15 0.1 0 Oct 120 100 80 60 0 Oct Nov Nov Dec Dec Jan Jan Feb Feb Mar 66 Apr Mar Apr Apr Snow Melted Mar Mean Tair>0 C Feb o Jan Snow Melted Mean Tair>0 C o Mean Tair~0 o C Mean Tair~0 o C Mean Tair~0 o C FrostDepth>80cm Tsoil2cm<0 o C -18 -20 10.2GHz 13.3GHz 16.7GHz 16.7GHz - 10.2GHz May May May 40 20 May Jun 200 -160 -120 Jun SoilW-2cm 0.05 Jun Dmax (mm) 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 Jun 2 6 GHz - -16 4 VV 16.7 0 8 VV 10.2 0 Snow Melted Mean Tair>0 o C Mean Tair~0 o C Mean Tair~0 o C Mean Tair~0 o C FrostDepth>80cm Tsoil2cm<0 o C 12 10 GHz (dB) IOP2 2010-2011 50 0 -50 Snow Depth (cm) Apr Snow Melted 0.2 Dec Mar Mean Tair>0 C 0.25 Nov Mean Tair~0 o C -200 Oct 0.3 Feb o -80 Jan Mean Tair~0 o C 0 Mean Tair~0 o C Mean Tair~0 o C -40 Dec Mean Tair~0 o C Mean Tair~0 o C 40 FrostDepth>80cm 80 Nov FrostDepth>80cm 120 Tsoil2cm<0 o C -22 Oct Tsoil2cm<0o C 160 Snowcover Onset -14 Snowcover Onset VV (dB) 0 -12 Snowcover Onset Soil Frost Depth (cm) -10 Snowcover Onset Soil Water Content (m 3 /m 3 ) -8 14 0 Jul 100 -100 Jul 1 0.8 0.6 0.4 0.2 0 Jul 120 100 80 60 40 20 0 Jul Figure 17. Measured backscattering coefficient, snow depth & soil frost depth, soil moisture content and snow grain size (Dmax) (from top to bottom) at Sodankyla during IOP2. Snow Depth (cm) -40 0 -80 0 Oct Nov Nov 120 100 80 60 40 Nov Dec Dec Dec Jan Jan Jan Feb Feb Mar Mar Feb Mar 67 Apr Apr Apr Snow Melted Mean Tair>0 o C Tsoil80cm<0 o C SoilW Stopped Decreasing Tsoil40cm<0 o C 10.2GHz 13.3GHz 16.7GHz 16.7GHz - 10.2GHz 2 May 20 May Jun 200 -160 -120 May May Jun Jun Dmax (mm) 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 Jun 6 GHz - -16 4 VV 16.7 0 Snow Melted Mean Tair>0 o C Tsoil80cm<0 o C SoilW Stopped Decreasing Tsoil40cm<0 o C Tsoil10cm<0 o C Tsoil2cm<0 o C Tsoil5cm<0 o oC Tsoil20cm<0 C 8 VV 10.2 0 12 10 GHz (dB) IOP3 2011-2012 50 0 -50 Snow Depth (cm) Apr Snow Melted 40 Mar Mean Tair>0 o C 80 Feb Tsoil80cm<0 o C 120 Jan SoilW Stopped Decreasing 160 Dec Tsoil40cm<0 o C -200 Oct Oct Nov Tsoil10cm<0 o C -22 Oct Tsoil10cm<0 o C -14 Snowcover Onset -12 Tsoil2cm<0 o C Tsoil5cm<0 o oC Tsoil20cm<0 C VV (dB) 0 -10 Tsoil2cm<0 o C Tsoil5cm<0 o oC Tsoil20cm<0 C -20 Snowcover Onset -18 Snowcover Onset Soil Frost Depth (cm) -8 14 0 Jul 100 Jul Jul -100 120 100 80 60 40 20 0 Jul Figure 18. Measured backscattering coefficient, snow depth & soil frost depth, soil moisture content and snow grain size (Dmax) (from top to bottom) at Sodankyla during IOP3. IOP4 2012-2013 Snow Melted 10 10.2GHz 13.3GHz 16.7GHz 16.7GHz - -20 Apr May 2 10.2GHz Jun -50 -100 Apr SoilW Stopped Decreasing May 40 20 0 Oct Nov Dec Jan Feb Mar Apr Jun Snow Melted Mar Tair>0o C Feb Mean Jan Tsoil80cm<0 o C Dec Tsoil2cm<0 o C Tsoil5cm<0 oCo Tsoil10cm<0 C Snowcover Onset Snow Depth (cm) 60 Nov Tsoil20cm<0 o C Tsoil40cm<0 o C -200 Oct 80 Snow Depth (cm) Snow Melted 0 -160 100 0 Jul 50 -120 120 4 100 Mean Tair>0 o C -80 Mar Tsoil80cm<0 o C 0 -40 Feb SoilW Stopped Decreasing 40 Jan Tsoil20cm<0 o C Tsoil40cm<0 o C 80 Dec Tsoil2cm<0 o C Tsoil5cm<0 o C Tsoil10cm<0 o C 120 Snowcover Onset Soil Frost Depth (cm) 160 Nov GHz - -18 6 VV 16.7 0 8 -16 -22 Oct 200 GHz (dB) 12 VV 10.2 0 Mean Tair>0o C 14 Tsoil80cm<0 o C SoilW Stopped Decreasing -14 Tsoil20cm<0 o C Tsoil40cm<0 o C VV (dB) 0 -12 Tsoil2cm<0 oo C Tsoil5cm<0 C Tsoil10cm<0 o C -10 Snowcover Onset -8 May Dmax (mm) 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 Jun Jul 120 100 80 60 40 20 0 Jul Figure 19. Measured backscattering coefficient, snow depth & soil frost depth, soil moisture content and snow grain size (Dmax) (from top to bottom) at Sodankyla during IOP4. Therefore, if the simulations of snowpits over partly-frozen soil and the snowpits influenced by snowmelt or instrument reinstallation are excluded, the mean bias and the 68 rms error of the MEMLS3&a simulations can be largely reduced. Table 7 shows the statistics for simulations between January 11, 2010 and May 19, 2010 for IOP1, and between December 26, 2010 and February 28, 2011 for IOP2, and between February 12, 2012 and March 1, 2012 for IOP3, and between January 19, 2013 and February 26, 2013 for IOP4. The starting date of the statistics is the time when the frost soil depth exceeded 80 cm and the soil moisture at 80-cm stopped decreasing; The ending date is when the mean air temperature first raised to about 0 ˚C with a frequent fluctuation in the SnowScat measurements due to snowmelt, or when the instrument was reinstalled. Time IOP1 IOP2 IOP3 IOP4 Freq. (Pol.) Mean bias (dB) VV VH Rms error (dB) HH VV VH HH Number 10.2 -0.51 0.46 2.97 0.64 0.64 2.99 13.3 0.07 0.77 3.00 0.64 1.00 3.02 16.7 -0.82 -0.66 0.75 1.09 0.94 1.01 10.2 0.37 0.92 3.29 0.56 1.02 3.31 13.3 0.03 0.75 2.98 0.52 0.93 3.01 16.7 -2.26 -0.80 -0.07 2.33 0.99 0.38 10.2 -1.42 -0.81 1.84 1.42 0.81 1.84 13.3 -0.20 0.14 4.06 0.20 0.14 4.06 16.7 0.06 0.98 3.20 0.06 0.98 3.20 10.2 -0.56 0.59 2.78 1.17 1.12 2.88 13.3 -1.14 -0.03 1.81 1.67 1.29 2.03 16.7 -0.67 -0.24 1.85 1.45 1.34 2.13 12 7 1 5 Table 7. Mean bias and rms error of MEMLS3&a simulations for snowpits not influenced by part-frozen soil conditions. As shown in Table 7, the rms error of VV- and VH-pol. backscattering coefficient simulations after partly-frozen soil and snowmelt are excluded reduce about 50% in most cases. The smallest rms error reached 0.52 dB, and is close to the +/- 0.5-dB error of the 69 SnowScat instrument. Because no pec optimization was applied, the average rms error is still about 1 dB for all the VV- and VH-pol. simulations. It may have some influences in the SWE retrieval results using MEMLS3&a for the SnowScat measurements, as will be studied in Chapter 4. 70 Chapter 3: SWE Retrieval Using Passive Microwave Measurements In this chapter, the SWE retrieval algorithm using passive microwave TB measurements is presented. The MEMLS-IBA was most accurate among the models and volume scattering coefficient parameterizations considered in Chapter 2.3, and therefore is used as the observation model here. The goal is to build a physically-based SWE retrieval algorithm allowing for multi-layer snowpack and prior information on the snow properties. The Bayesian Algorithm for SWE Estimation with Passive Microwave measurements (BASE-PM) is used as the retrieval algorithm. It is based on the Markov Chain Monte Carlo (MCMC) method (Gelman et al., 1995). Two kinds of priors were prepared and tested. The first set of priors is from globally available datasets, hereafter referred as “generic priors” because they represent estimates of snow parameters and their uncertainty that are available globally. The second set of priors is from historical local snowpit measurements in the area, referred as “local priors”. The specific research questions are: (1) Can BASE-PM using the generic prior provide SWE estimates that meet the Integrated Global Observing Strategy (IGOS, 2007) requirement for shallow snow, which is 30 mm root mean squared (RMS) error? (2) How does the SWE prior influence SWE estimate accuracy? 71 (3) How does the prior for the stratigraphy information (layered density, temperature and grain size) influence the accuracy of the estimated SWE? Section 3.1: Estimated Parameters for Passive SWE Retrieval The brightness temperature (TB) observed at the snow surface is basically determined by two sets of snowpit parameters: the first set is the snow parameters, including snow depth (thickness), micro-structure parameter (Dmax, pec or others), snow density, snow temperature and snow wetness; the second set is the soil parameters, including soil moisture, soil temperature, and other soil texture-related parameters (including soil bulk density, and sand, silt and clay contents). Therefore, to retrieve SWE, not only SWE itself, but other parameters need to be estimated at the same time. Based on the large ratio of the wavelength relative to the geometry of the snow microstructure, volume scattering effect dominates the microwave response for dry snowpacks compared to absorption effect (Wiesmann and Matzler, 1999; Liang et al., 2008; Lemmetyinen et al., 2010; Picard et al., 2013). Volume scattering strongly decreases the TB observed at high frequencies, with an extent determined by the number of scatters (ice particles in snow) and the dimension of the scatters (Ulaby et al., 2014). Therefore, among the snow parameters, the estimation of snow depth and snow micro-structure parameter are of critical importance. Other snow parameters need to be estimated are snow density and temperature. Snow density influences the diffraction angle of microwave when it passes through the air-snow boundary; it also slight influences the volume scattering coefficient through the snow dielectric constants. Snow physical 72 temperature determines the internal thermal source of snow. In general, their influences to TB are much smaller than the snow thickness and snow micro-structure parameter (Durand et al., 2008), and thus their estimation may partly rely on the prior information. Snow wetness is not considered for dry snow. The other issue is, do we have to consider multi-layer snow effects. Typically, for homogenous snow, TB at the snow surface decreases with increasing SWE and gradually becomes saturated for thick snow (Pan et al., 2016). However, Matzler et al.(1982) and Derksen et al. (2010) observed an increasing TB trend after the saturation point for natural snowpacks. One hypothesis is, it can be related to the vertically inhomogeneous snow properties inside natural snowpacks. Commonly, older snow with depth hoar and refrozen crusts is located beneath new snow with smaller grain size (Lemmetyinen et al., 2016). In this case, microwave radiation with different penetration depths may observe different snow properties in the snowpacks. Therefore, when BASE-PM is used with a multi-layer snow RT model, 1- up to 6-layer estimation for snow properties were tested. The increased number of estimated parameters challenges the ability of limited microwave measurements to estimate all of them. In this case, the prior information is required, as mentioned in Durand & Liu (2012). Soil properties provide the background radiation prior to the snow attenuation. In this research, soil is simply considered as the half-space medium with rough surface. The soil-emitted TB is determined by the soil temperature and the soil surface emissivity. Here, the rough bare soil reflectivity model (Wegmuller and Matzler, 1999) is used to calculate soil emissivity (i.e., 1- reflectivity) from the soil surface roughness and soil 73 dielectric constants. Hereafter, we call it the WM model. The soil roughness parameter used in WM is the rms-height of the soil surface. The Dobson semi-empirical soil dielectric constant model (Dobson et al., 1985) is used to calculate the soil dielectric constants. The Dobson model requires soil moisture, soil temperature, soil bulk density, and soil sand, silt, clay contents to calculate the soil dielectric constants. Here, a simplification is made to estimate only the soil roughness, soil temperature and soil moisture in BASE-PM. The soil bulk density is fixed as 1.5 g/cm3; soil texture is fixed as 30.63% sand, 55.89% silt and 13.48% clay, which is a medium texture for experiments in Hallikainen et al. (1985). The original Dobson model does not treat frozen soil. Thus, the estimated soil moisture must be considered to be an effective value if soil temperature is below the freezing point (Mironov et al., 2010). Section 3.2: The BASE-PM Retrieval Algorithm Theory and Implementation The Markov Chain Monte Carlo (MCMC) method is based on Bayesian statistical theory, and is a numerical realization of Bayes’ law. When MCMC is applied to estimate SWE, it aims to find a combination of estimated snowpack variables that gives the highest possible posterior probability conditioned on the TB measurements and the prior information. In other words, if a vector containing all the estimated variables is θ and a vector containing the microwave measurements is y, the goal of SWE retrieval is to find the highest p(θ|y). According to Bayes’ law, p(θ|y) can be calculated as (Gelman et al., 1995): | = ! ! ! !|! (3.1) ! ! 74 In this equation, = ! | , which is the sum of TB over all possible θ. p(y) is invariant for a particular snowpit, and thus its calculation can be avoided, as described below. p(θ) is the prior probability of the snow and soil parameters for a particular snowpit. p(y|θ) is the probability of the observed y given θ: it is essentially just a statistical model of the observation error. Assumed the TB observation error follows a multi-variate normal distribution, p(y|θ) can be calculated as (Durand & Liu, 2012): | = 2 ! !! ! Σ! ! ! ! ! − ! − ! !! Σ! − (3.2) where Ny is the number of TB measurements in y, Σ! is the TB error covariance matrix (2 K diagonal matrix), and M(θ) is the observation model, which is MEMLS-IBA with WM soil reflectivity model and Dobson soil dielectric constant model. Using equation (3.2), the snow microwave emission physics is involved in retrieval. The MCMC algorithm computes p(θ|y) by starting from an initial combination of snowpack variables, and using a “likelihood ratio” to guide their estimates into the region where they best represent the observed TB via the observation model simulation. As shown in Figure 20, at iteration i, MCMC “randomly walked” to a new sets of θi,candidate from θi-1. A symmetric jump function is used to conduct this random walk, as J(θi) = θi-1 + dθi, where, dθi, is a normally-distributed random number centered at zero with a standard deviation referred to the jump size. The jump size is adaptively tuned to provide the optimal acceptance rate (Gelman et al., 1995). Then, p(θi-1| y) and p(θi,candidate|y) are compared as: = ! !!,!"#$%$"&' |! ! !!!! |! = ! !!,!"#$%$"&' ! !|!!,!"#$%$"&' ! !!!! ! !|!!!! 75 (3.3) Note that p(y) has been eliminated. r is called the likelihood ratio (Gelman et al., 1995). If r is rc-percent probably larger than 1, where rc is a probabilistic value, θi,candidate is more likely better than θi-1 and will be accepted as θi. Otherwise, θi=θi-1, and the algorithm will try another jumping direction in the next iteration. rc instead of 1 is used to prevent local optimization in the Markov Chain realization. ith iteration i=i+1 Go to next iteration θi,candidate = θi-1 + dθi θ =θ i p(θ i,candidate |y) p(θi-1 |y) r r >= rc i, candidate θi=θi-1 r < rc Figure 20. A schematic view of the flow chart of the MCMC algorithm at the ith iteration. The iteration i starts from the up left corner, where a candidate combination of snow/soil parameters, θi,candidate is produced using jump functions from θi-1. θ includes multi-layer snow thickness, temperature, density, exponential correlation length, and soil temperature, soil volumetric moisture and soil surface roughness. Then, the posterior probability p(θi,candidate |y) and p(θi-1 |y) were calculated and compared using equation (3.3). The algorithm will determine final estimate at i-th iteration (θi) according to r with relevant to a probabilistic value, rc. After a large number of iterations, a Markov chain will be built. The chain numerically sampled the posterior snow and soil parameters that can reproduce the microwave TB measurements, except for the first several iterations that are still influence by the initial values. These iterations are called the burn-in period, and were removed from the calculation of final MCMC estimates. Taking SWE as an example, the MCMC-estimated SWE is calculated as: 76 !"!" = ! ! !"!#$ !!!"#$!!" !!"!#$ !!!!"#$!!" !! ! (3.4) where, SWEi is SWE at the i-th iteration. itotal is the total length of the MCMC chain. iburnin is the length of the burn-in period. The uncertainty of MCMC-estimated SWE is calculated as: !"#,!"!" = ! !!"!#$ !!!"#$!!" !!"!#$ !!!!"#$!!" !!(! − !"!" )! (3.5) The MCMC-estimates of other variables can be calculated using similar methods. Section 3.3: Priors for SWE Retrieval 1) Local Priors The local priors are computed from statistics of snowpit measurements. For SWE, the local prior was computed from the monthly-average SWE of all snowpits at each site. The local priors for density, temperature and grain size metrics were computed from layered snow property measurements. It should be noted that, MCMC algorithm was run for 1-, 2-, up to 6-layer snow estimates before it can determine which number of layers (nlyr) gives highest posterior probability. For prior calculation, if the actual number of snow layers is different from 1, 2,.. or 6 layers, it was first relayered to the required number of layers using the method in Jordan (1991). Basically, the same datasets used in RT model comparison in Section 2.3 were used for SWE retrieval in this chapter. To prevent the local prior from too accurate, some originally-utilized snowpits were reserved for prior calculation, or some other snowpit 77 measurements at the nearby location or from other years were added. See details as follows: (1) At Sodankyla, SWE retrieval was conducted on 36 snowpits measured in IOP 1 and 2 (2009/2010, 2010/2011 winters), using TB measurements at 10.65, 18.7, 36.5 and 90 GHz, vertical polarization at 50° observation angle (Lemmetyinen et al., 2016). Snowpits in IOP 3 and 4, and other snowpits in IOP 1 and 2 with missing data were reserved for local prior calculation. This is partly because in IOP 3 and 4, 23-GHz radiometer was operated instead of 90-GHz. (2) At Churchill, 6 snowpits measured in forest-clearings by Derksen et al. (2012) were used for SWE retrieval (the other 2 has no TB measurements at 6.925 or 89 GHz). More snowpits observed in Churchill area during the Canadian Cold Regions Hydrology HighResolution Observatory (CoReH2O) Snow and Ice Experiment and the follow-on experiments (Montpetit et al., 2013) were used for prior calculation. These snowpits include some deep taiga snow, tundra snow and snow over dry/wet fen and lake ice. (3) At Colorado during CPLX experiment, 6 snowpits observed at LSOS during IOP 3 on the week of February 19-26, 2003 (Cline et al., 2003) were used for SWE retrieval. The snow type is alpine snow. The local priors were calculated from 73 ISA snowpits from Fraser MesoCell Study Area (MSA) where LSOS located. Table 8 lists the mean and the standard deviation of SWE and mass-weighted average density, temperature and exponential correlation length (pec) for all snowpits at different months. pec was converted from the geometric snow grain size (Dmax) using the same methods shown in equation (2.33). It can be seen from Table 1 that, the standard 78 deviation of SWE varies from 17% to 41% of the mean SWE at different months and sites. However, note that sample sizes are rather small, in some cases. Therefore, for experimental purpose, the standard deviation of the local SWE prior is determined as 40% of the measured mean SWE. In other words, a 40% relative accuracy for local SWE prior was assumed for all cases. Site, Month CLPX Churchil l Sodanky lä Feb Jan Feb Mar Apr Nov Dec Jan Feb Mar SWE (mm) Mean 185.08 112.42 137.19 179.14 149.99 54.71 75.82 110.29 150.13 173.83 Std 75.67 22.64 36.96 70.93 51.90 13.27 21.99 25.31 30.91 29.12 Density (kg/m3) Mean 222.85 239.82 279.91 296.53 321.64 222.03 219.98 218.59 234.30 227.64 Tsnow (°C) Std 35.71 29.78 57.04 31.18 50.80 29.66 31.88 30.72 36.12 18.95 Mean -4.12 -8.86 -11.04 -9.68 -3.28 -3.39 -3.06 -5.89 -8.04 -6.23 pec (mm) Std 1.38 2.37 4.62 5.91 2.07 3.11 3.33 2.21 4.44 3.10 Mean 0.174 0.203 0.236 0.214 0.206 0.175 0.205 0.236 0.238 0.158 N Std 0.041 0.015 0.041 0.053 0.052 0.049 0.042 0.041 0.032 0.044 73 5 17 16 7 8 27 21 18 33 Table 8. Statistics of SWE and the mass-weighted average density, temperature and micro-structure parameters of the measured snowpits. 2) Generic Priors Unlike the local prior, the generic prior is from a set of globally available datasets. The generic prior for SWE is from the global Variable Infiltration Capacity (VIC) Macroscale Hydrologic Model predictions (Nijseen et al., 2001). The VIC model uses the timeseries meteorological measurements to estimate the evolution of snowpack, including snow accumulation from snowfalls, snow compaction and grain growth, and snow ablation due to snowmelt, driven by surface mass and energy balances. The spatial resolution is 2-degree, and the time 79 range is from 1980 to 1993 (http://vic.readthedocs.org/en/master/). Monthly-averaged 2-degree VIC SWE was used for prior estimate at Sodankylä and Churchill. At LSOS, however, it was found that a 2degree resolution is too coarse to capture the snow accumulation in mountain. The 2degree VIC model gives a SWE of 6.6±6.2 mm in February, but the mean SWE from snowpit measurements is 185.08 mm (see Table 1). In this case, a 0.125-degree VIC SWE product driven by the North American Land Data Assimilation System (NLDAS), Phase 2 is chosen for the CPLX-LSOS site. The data is downloaded from NLDAS-2 Monthly Climatology datasets (Mocko, 2014). Table 9 lists the monthly-average SWE for different months. Comparison between Table 8 and Table 9 shows that, the SWE predicted by VIC model is underestimated by about 30% to 58% compared to the locally-measured SWE. Therefore, the standard deviation of the generic SWE prior is set as 100% of the VIC SWE. The ratio of standard deviation to the mean is higher for the generic prior, reflecting the lack of precision currently available for global SWE estimates. 80 Site, Month VIC SWE (mm) Multi-year average (mm) CLPX Churchill Sodankylä Feb Jan Feb Mar Apr Nov Dec Jan Feb Mar 90.52 77.20 95.79 108.37 104.59 23.01 41.17 65.57 89.02 115.45 Standard deviation between different years (mm) 38.76 25.47 24.41 25.69 35.78 17.93 34.91 52.13 58.49 58.98 Generic density prior (kg/m3) Mean Std 335 217 86 56 Generic Tsnow prior (°C) Generic pec prior (mm) Mean Mean -6.5 -10.5 Std 6.5 10.5 0.18 Std 0.18 Table 9. Monthly-averaged SWE from the VIC model predictions and the generic priors of snow parameters used for retrieval The generic priors for snow density and temperature are determined according to typical snow classification system (Sturm et al., 1995). At Sodankylä and Canada, the snowpits we used for retrieval are classified as taiga snow, and the snowpits at CLPX-LSOS are classified as alpine snow. Therefore, the priors are chosen strictly according to the snow type. Sturm et al. (2010) presents a bulk density of 217±56 kg/m3 for 1541 taiga snowpits, and a bulk density of 335±86 kg/m3 for 4623 alpine snowpits. These values were used directly as mean and standard deviation of snow density to calculate the lognormal distribution parameters (Pan et al., 2017). Sturm et al. (1995) shows the lower bound of the air temperature for taiga snow is -21 °C, and the upper bound of the soilsnow interface temperature is 0 °C. Therefore, we assumed the prior of snow temperature for taiga snow is -10.5±10.5 °C. The prior of snow temperature for alpine snow is 6.5±6.5 °C, as determined by the same method. 81 For the snow micro-structure information not provided by Sturm’s classes, the generic prior for pec was assumed to be 0.18±0.18 mm. It was determined according to the measurements published in Matzler (2002), where pec of fresh snow is about 0.03 mm, and pec of depth hoar is about 0.32 mm. As a log-normal distributed variable, a mean and standard deviation of 0.18 mm for pec correspond to a range of 0.05-0.29 mm if one standard deviation around the mean of log(pec) is considered. According to Durand et al. (2008), a 0.18-mm pec is approximately equivalent to a 1-mm geometric grain size (Dmax). Generic priors for snow parameters were the same for different snow layers. However, the surface snow density was constrained to be smaller and the surface snow temperature was constrained to be lower than the bottom layers during the MCMC iterations. 3) Other Priors As to the other priors, the prior for number of layers (nlyr) follows a Poisson distribution. The parameter λ of the Poisson distribution is 2 for taiga snow, and 5 for alpine snow. It is determined by the fact that usually alpine snow is thicker and thus has the potential to have more vertical variations. The prior for soil temperature was set to be the same as the prior used for the temperature of the bottom-most snow layer. The prior for soil moisture was fixed at 8±8%, and the prior for soil roughness was fixed as 1±1 cm, following lognormal distributions as the snow variables. 82 Sodankylä Churchill CLPX-LSOS (a) Measurement 0.8 0.7 109 kg/m , 0.09 mm 3 z (m) 0.6 165 kg/m3, 0.09 mm 0.5 224 kg/m3, 0.28 mm 0.4 (b) Local priors 205±22 kg/m3, 0.14±0.04 mm (c) Generic priors 217±56 kg/m3, 0.18±0.18 mm 199 kg/m , 0.28 mm 3 0.3 0.2 245 kg/m3, 0.28 mm 247±24 kg/m3, 0.31±0.03 mm 0.1 276 kg/m3, 0.28 mm 217±56 kg/m3, 0.18±0.18 mm 0 Figure 21. The in-situ measurements (a) and the priors (b)-(c) for a snowpit at Sodankyla observed on March 2nd, 2010 as an example. The local priors used the average values of the measured snow properties at Sodankyla in March, see (b). The generic priors used the VIC-model SWE predictions and the Sturm’s snow classes; the generic prior of pec is 0.18±0.18 mm, see (c). The priors for layer thickness were converted from SWE and density priors. On the up-right corner, it shows the geographic location of the sites used in this paper. 4) An Example of Priors All the estimated variables were assumed to follows log-normal distributions as in Durand & Liu (2012). The method to calculate lognormal distribution parameters and to convert SWE prior to snow layer thickness prior can be found in Pan et al. (2017). 83 Figure 21 shows an example of the local and generic priors for a snowpit observed at Sodankylä on March 2nd, 2010. Figure 21(a) shows the in-situ snow measurements. The measured snowpit has 6 layers. The density is 109.1 kg/m3 at the surface layer, and it increases to 275.75 kg/m3 at the bottom layer. The pec converted from measured Dmax is 0.09 mm for the top two layers, and is 0.28 mm for the other layers. The measured SWE is 165.9 mm. Figure 21 (b) shows the local priors using 2 layers. The local SWE prior is 173.84 mm, which is 4.7% larger than the true SWE. It was converted into slightly shallower snow depth, because the local prior of density is larger than the measured density of this snowpit. Figure 21 (c) shows the generic priors using 2 layers. The generic SWE prior is 115.45 mm, which is 30.4% smaller than the true SWE. It was converted in to a snow depth 36.7% shallower than the measured depth. Section 3.4: Passive Microwave Retrieval Results In this section, the MCMC Retrieval algorithm was tested over difference choices for SWE priors, stratigraphy priors (i.e., priors for snow temperature, density and pec) and the inputted TB observations. There are two types of SWE priors: local prior from measured mean SWE, and generic prior from VIC model predictions. There are two types of stratigraphy priors: local prior from the statistics of snowpit measurements, and generic prior from snow classes. Besides the actual TB measured by the ground-based radiometers, the measured snow and soil properties for each snowpit were also used to generate a synthetic TB calculated by the same observation model used for retrieval. A 2K random error is added to the synthetic TB, where 2 K is the observation error of the 84 instrument. The synthetic experiments are presented in order to help understand and explain the algorithm performance using the real TB. Retrieved SWE Using Real TB Measurements Figure 22 shows the performance of the SWE retrieval using real TB measurements with generic SWE and stratigraphy priors, where (a) is the BASE-PM estimated SWE versus the true SWE, and (b) is the BASE-PM estimated SWE uncertainty versus absolute SWE error. This is the result for all the 48 snowpits from Sodankylä, Churchill and CLPXLSOS. BASE-PM SWE performance statistics are shown in Table 10, row 1. For all 48 snowpits, the rms error of estimated SWE for all snowpits is 57 mm, and the mean bias is -9.3 mm. However, as shown in Figure 22(a), there are two snowpits from the Churchill site, which have significantly higher SWE estimate than the true SWE. A Grubbs’ test (Grubbs, 2007) applied on the absolute SWE error indicates these two snowpits are outliers with a significance level of 0.005. In Figure 22 (b), it shows the uncertainty estimates are also highest for these two snowpits; thus it is reflected from the retrieval algorithm that their SWE estimates are unreliable. After excluding the two outliers, the rms error of the remaining 46 snowpits is 42.7 mm, and the mean bias is -17.3 mm. They are improved compared to the rms error (52.7 mm) and the mean bias (-42.1 mm) of the generic SWE prior. In Figure 22(a), the SWE of all CLPX snowpits are underestimated. The reason for the underestimation is most likely due to the reduced sensitivity of passive microwave signal to SWE for thick snow. Another reason could be only TB at three frequencies (18.7, 36.5 and 89 GHz) was measured compared to the four frequencies 85 utilized at two other sites (10.65, 18.7, 36.5 and 90 GHz at Sodankyla and .9, 19, 37 and 89 GHz at Churchill). The frequency choice will be checked later in the synthetic experiments. If both the outliers and the CLPX snowpits are excluded, the rms error of BASE-PM is reduced to 30.8 mm, and the mean bias is reduced to -7.1 mm. In this case, the algorithm meets (or is within 0.8 mm) of the 30-mm Integrated Global Observing Strategy (IGOS) requirement for shallow snow (IGOS, 2007). It should be noted that the “shallow snow” defined in IGOS is <300 mm SWE, whereas the rms error of our algorithm summarized here is for SWE<170 mm. (a) (b) 400 300 250 200 BASE−PM SWE uncertainty (mm) sd−nov sd−dec sd−jan sd−feb sd−mar ch−jan ch−feb ch−mar ch−apr clpx−feb 350 SWEBASE−PM (mm) 300 150 100 50 0 0 100 200 SWEtrue (mm) 300 400 250 200 150 100 50 0 0 100 200 |SWE error| (mm) 300 Figure 22. The performance of the BASE-PM (Bayesian Algorithm for SWE EstimationPassive Microwave) for SWE retrieval using real TB measurements: (a) the BASE-PM SWE compared to the observed SWE at different sites and months, (b) the BASE-PM SWE uncertainty compared to the absolute SWE error. BASE-PM used the generic priors for both SWE and the snow stratigraphy feature. 86 In Figure 22 (b), it also shows the uncertainty estimates clustered around 25-50 mm, which is indeed the precision of the method, excluding the outliers and the CLPX snowpits. The uncertainty for the CLPX snowpits is underestimated; this may be due to underestimation of the prior uncertainty. The fact that the two outlier snowpits have the highest uncertainty suggests that the BASE-PM uncertainty estimate is able to flag poor retrievals. More results using other prior combinations can be found in Figure 23, with statistics of their performance shown in Table 10, row 2-4. As mentioned before, the SWE rms error using generic SWE and stratigraphy priors was 42.7 mm without outliers. In Figure 23, the RMS error using generic SWE prior and local stratigraphy prior for the same snowpits is 33.3 mm. The RMS error using local SWE prior and generic stratigraphy prior is 29.0 mm. And, finally, the RMS error using both local priors is 28.3 mm. The statistics without CLPX snowpits is also available. Compared to when the generic priors were used, the rms error decreased 22% (18% without the CLPX snowpits) when a better stratigraphy prior was introduced. The rms error decreased 32% (22% without the CLPX snowpits) when a better SWE prior was introduced. The rms error decreased 34% (26% without the CLPX snowpits) when both a better stratigraphy prior and a better SWE were introduced. For all cases, the mean bias of SWE estimates is smaller than the local SWE prior. The rms error of SWE estimation improves over the rms of local SWE prior when the outlier snowpits and the CLPX snowpits are excluded. Thus, prior information on either SWE or stratigraphy improves the BASE-PM SWE estimation performance. 87 SWEMCMC(mm) 300 sd−nov sd−dec sd−jan sd−feb sd−mar ch−jan ch−feb ch−mar ch−apr clpx−feb local SWE prior generic Strati. prior 200 100 0 0 100 200 300 SWEtrue (mm) MCMC SWE uncertainty (mm) SWEMCMC(mm) 300 generic SWE prior local Strati. prior 200 100 0 0 100 200 300 SWEtrue (mm) MCMC SWE uncertainty (mm) 300 400 (e) 400 SWEMCMC(mm) 400 (c) 500 400 local SWE prior local Strati. prior 200 100 0 0 MCMC SWE uncertainty (mm) (a) 400 100 SWE 200 300 (mm) true 400 300 (b) 250 200 150 100 50 0 0 300 100 200 |SWE error| (mm) 300 (d) 250 200 150 100 50 0 0 300 100 200 300 |SWE error| (mm) 400 (f) 250 200 150 100 50 0 0 100 200 |SWE error| (mm) 300 Figure 23. The MCMC-estimated SWE compared to the observed SWE using observed TB from radiometers, when (a) local SWE prior but generic stratigraphy (Strati.) prior were used; (c) generic SWE prior but local stratigraphy prior were used; (e) local SWE prior and local stratigraphy prior were used. In (b), (d) and (f), it shows the MCMC SWE uncertainty versus the absolute SWE error, corresponding to (a), (c) and (e), respectively. 88 Experiment # Retrieved SWE Prior SWE 1 2 3 4 All snowpits SWE prior Stratigraphy prior Generic Generic Generic Local Local Generic Local Local Generic SWE prior Local SWE prior RMS error (mm) 57.0 64.2 34.1 36.0 52.3 25.5 Mean bias (mm) -9.3 3.7 -0.04 2.0 -42.0 11.7 Excluding two pits at Churchill RMS error (mm) 42.7 33.3 29.0 28.3 52.7 24.9 Mean bias (mm) -17.3 -7.6 -3.7 -2.5 -42.1 11.5 Excluding two pits at Churchill, and CLPX snowpits RMS Mean error bias (mm) (mm) 30.8 -7.1 25.2 0.8 24.1 3.0 22.8 4.6 36.9 -32.2 25.4 15.3 Table 10. Error of MCMC-retrieved SWE using observed TB Retrieved SWE Using Synthetic TB Figure 24 shows the retrieved SWE using the synthetic TB, and the statistics is shown in Table 11. First, In contrast with the results using real TB, the SWE outliers for these two Churchill snowpits are no longer found. The SWE RMS error and bias become similar with or without these two snowpits. Therefore, it shows the outliers may be related to some special features in the snowpack or the soil, or a limitation of the microwave model, and thus was not reflected in the synthetic TB. The underestimations for the CLPX snowpits were reproduced by the synthetic test. Also, it shows that the underestimation is not related to the frequency choice, because the results using 4-frequency synthetic TB with an additional frequency as 10.65 GHz (in black stars) are not significantly better than the results using 3-frequency synthetic TB (in black diamonds). In the synthetic test, the rms SWE error is also smaller when the more accurate local SWE or stratigraphy prior than the generic priors is used. It supports the conclusions found for tests using real 89 TB: using more accurate priors, where available, can improve the accuracy of SWE retrieval. (a) (b) 300 400 sd−nov sd−dec sd−jan sd−feb sd−mar ch−jan ch−feb ch−mar ch−apr clpx−feb clpx−feb2 generic SWE prior generic Strati. prior 200 SWEMCMC(mm) SWEMCMC(mm) 400 100 0 0 100 200 300 SWEtrue (mm) 300 200 100 0 0 400 local SWE prior generic Strati. prior 100 (c) 300 400 local SWE prior local Strati. prior generic SWE prior local Strati. prior SWEMCMC(mm) SWEMCMC(mm) 400 400 200 100 0 0 300 (d) 400 300 200 SWEtrue (mm) 100 200 300 SWEtrue (mm) 400 300 200 100 0 0 100 200 SWEtrue (mm) Figure 24. The MCMC-estimated SWE compared to the observed SWE using synthetic TB, when (a) generic SWE prior and generic stratigraphy (Strati.) prior were used; (b) local SWE prior but generic stratigraphy prior were used; (c) generic SWE prior but local stratigraphy prior were used; (d) local SWE prior and local stratigraphy prior were used. 90 Retrieved SWE Experiment # SWE prior choice Stratigraphy prior choice 5 6 7 8 Generic Generic Local Local Generic Local Generic Local All snowpits RMS error (mm) 47.6 33.9 28.4 26.4 Bias (mm) -29.5 -12.0 -1.9 -2.9 Excluding two pits at Churchill RMS error (mm) 46.7 34.5 28.6 25.8 Bias (mm) -28.2 -12.8 -1.2 -2.3 Excluding two pits at Churchill, and CLPX snowpits RMS Bias error (mm) (mm) 39.9 -21.1 24.5 -4.7 24.8 5.3 18.8 4.9 Table 11. Error of MCMC-retrieved SWE using synthetic TB Markov Chains of Estimated Variables Figure 25 to Figure 28 explain more details of the BASE-PM algorithm when it was applied on the example snowpit mentioned in Figure 21. Figure 25 shows the Markov chains of the snow layer thickness, density, pec and temperature when the generic SWE prior and stratigraphy prior were used. The chains randomly jumped within a large range for each variable after adjusted during the burn-in period. For example, in the first 2000 iterations of the burn-in period, the pec of the bottom layer increases from the prior value as 0.18 mm to above 0.3 mm. Then, it randomly fluctuates between 0.25 to 0.5 mm with a stable standard deviation. Using the values of variables at each iteration, a corresponding TB chain can be calculated, as shown in Figure 26. The RMS error of the TB chain compared to the observed TB is 2.19 K, 4.17 K, 2.07 K and 2.17 K at 10.65, 18.7, 36.5 and 89 GHz. The RMS error at 18.7 GHz is larger because the mean bias (-3.98 K) at this frequency is also larger than other frequencies (1.38 K, 0.76 K and -0.79 K at 10.65, 36.5 and 89 GHz, respectively). This may be due to some physical emission characteristics at 18.7 GHz that are not well91 represented by the observation model. Note that many other snowpits have a good (RMS error smaller than 2.0 K) TB chain at 18.7 GHz, too. The standard deviation of the TB chains all fall below 2 K for the four frequencies (1.71 K, 1.27 K, 1.93 K and 2.02 K at 10.65, 18.7, 36.5 and 89 GHz, respectively). This indicates that many combinations of snow/soil properties could produce a similar TB signal. However, if a SWE chain is produced using the layered snow thickness and density from these combinations, it could be found in Figure 27 that, the SWE chain basically follows a log-normal distribution. The uni-modal shape of the SWE chain indicates that there are some SWE values which have a higher probability than other values. This is because the MCMC algorithm sampled the combinations of snow variables conditioned on the TB observations. For the snowpit example here, the retrieved SWE is 144.3 mm with an estimated uncertainty of 42.9 mm. The true SWE is 165.9 mm, and thus the absolute error is 21.6 mm. The retrieved SWE is more accurate than the prior SWE estimate (115.45 mm). 92 (a) (b) 1.4 500 BotLyr SurfLyr 1 0.8 0.6 0.4 0.2 0 0 BotLyr SurfLyr 450 Layer density (kg/m3) Layer thickness (m) 1.2 400 350 300 250 200 150 100 0.5 1 1.5 Iteration # 2 2.5 50 0 3 0.5 1 4 x 10 (c) 2 2.5 3 4 x 10 (d) 275 Layer temperature (oC) 0.5 0.4 Layer pec (mm) 1.5 Iteration # 0.3 0.2 BotLyr SurfLyr 0.1 270 265 260 255 250 BotLyr 245 SurfLyr 0 0 0.5 1 1.5 Iteration # 2 2.5 240 0 3 4 x 10 0.5 1 1.5 Iteration # 2 2.5 3 x 10 4 Figure 25. The BASE-PM Markov chains of the snowpit observed at Sodankylä on March 2nd, 2010 for (a) snow layer thickness, (b) snow density, (c) pec, and (d) temperature, when generic SWE and stratigraphy priors were used. BASE-PM chose 2 layers, where “BotLyr” is for the chain of the bottom layer, and “SurfLyr” is for the chain of the surface layer. 93 280 260 TB (K) 240 220 200 180 10.65 GHz 18.7 GHz 36.5 GHz 89 GHz 160 140 0 0.5 1 1.5 Iteration # 2 2.5 3 x 10 4 Figure 26. Simulated TB from the estimated snow and soil variables in Figure 19, where the crosses are the measured TB at four frequencies. 2500 σSWE=42.9 mm bias=21.7 mm Frequency 2000 1500 1000 500 0 50 100 150 200 250 300 350 Snow water equivalent (mm) 400 450 Figure 27. Histogram of the SWE for the example Sodankylä snowpit, where the blue dash line indicates the BASE-PM SWE, and the back dash line indicates the measured SWE. 94 Figure 28 shows the histograms of the snow layer thickness, density, temperature and pec. The posterior distributions of these variables are also uni-modal. Table 10 summarizes the retrieved values for all snow and soil properties. The “true” 2-layer thickness, density, temperature and pec shown here were derived from the measured snow properties using the relayering approach mentioned before. The retrieved snow depth is 0.69 m, which is 17.9% shallower than the true snow depth, but the layer thicknesses are quite different from the relayered values. Thus, snow property estimates for individual layers are generally less accurate than bulk estimates. For example, it is also very interesting to find that the overall SWE accuracy is better than the accuracy of layer thickness or the density for specific layers. The retrieved layer thickness is 50.0% lower at the bottom layer and 25% higher than the surface layer; and, the retrieved density is 0.5% higher for the bottom layer and 37.7% higher for the surface layer. On the contrary, the retrieved SWE is only 13.1% less than the true SWE. However, the algorithm successfully estimated a lower pec at the surface layer than the bottom layer, when 0.18 mm pec prior was used for both layers. The reason for the less accurate snow properties at specific layers is most likely due to the lack of enough prior information to constrain the estimates. Of all the snow variables, the sensitivity of passive microwave TB to density and temperature is small. To match the TB measurements, the layer thickness and pec are most important, and their estimates can largely compensate for the effects in density and temperature estimates. However, there are still many combinations of layered thickness and pec that can give similar TB signals. For example, as shown in Fig. 25(a) and (c), the thickness of the bottom layer shows a slight negative correlation with the pec of the bottom layer. The 95 relatively large pec at the bottom layer in Fig. 25(c) corresponds to the relatively smaller thickness of the bottom layer. These two parameters compensate for each other. This highlights the importance of the histogram of posterior SWE; the histogram provides the uncertainty embedded in the uncertainty of all snow parameters that can influence the SWE estimation. (a) (b) 3000 2500 BotLyr SurfLyr 2000 2000 Frequency Frequency 2500 BotLyr SurfLyr 1500 1000 1500 1000 500 500 0 0 0.5 1 0 0 1.5 Layer thickness (m) 100 (c) 3 400 500 (d) 3000 BotLyr SurfLyr 5000 BotLyr SurfLyr 2500 4000 Frequency Frequency 300 Layer density (kg/m ) 6000 3000 2000 1000 0 0 200 2000 1500 1000 500 0.1 0.2 0.3 Layer pec (mm) 0.4 0 0 0.5 10 20 30 274 − Layer temperature (K) 40 Figure 28. Histogram of the chains in Figure 25 after the burn-in period for: (a) layer thickness, (b) density (b), (c) pec, and (d) 274 K minus snow temperature in K. In (d), the variable estimated related to temperature is 274 K minus the temperature in K, because the assumption of a log-normal distribution requires positive values. 96 Estimated Variables Snow layer thickness (m) MCMC estimates True values Bottom Layer 0.26 0.52 Surface Layer 0.40 0.32 Bottom Layer 242.85 241.57 Surface Layer 189.56 137.69 Bottom Layer 0.320 0.28 Surface Layer 0.051 0.09 Bottom Layer -4.4 -6.56 Surface Layer -8.9 -8.42 Snow depth (m) 0.69 0.84 Snow water equivalent (mm) 144.2 165.9 Soil temperature ((˚C) -4.21 -4.27 Soil moisture (%) 7.01 4.64 Soil roughness (cm) 0.7 0.1 Snow density (kg/m3) pec (mm) Snow temperature (˚C) Table 12. The estimated snow and soil parameters for the example Sodankyla snowpit observed on March 2, 2010 when generic SWE and stratigraphy priors were used. The “true” values of the snow variables were relayered from the original 6-layer measurement. For soil parameters, the soil temperature was estimated with a relative bias of -1.4%. The estimated soil roughness and moisture are both different with the true values. It was usually found these two parameters compensated with each other to produce a similar soil-emitted TB at low frequencies. This is because only TB at vertical polarization is not enough to determined two parameters, whereas TB at horizontal polarization is so strongly influenced by ice layers in snow that can not be used for SWE retrieval. Use of the TB measurement before snow season, if available, may help to determine soil roughness, but it requires the bare soil without vegetation; this will be explored in future work. 97 Section 3.5: Discussions In this section, the reasons for the two outlier snowpits at Churchill and the underestimated CLPX snowpits are studied by checking the response of multi-frequency TB to SWE from observations and simulation experiments. The MCMC retrieval algorithm was run for 1 up to 6 snow layers, however, for most snowpits, the MCMC algorithm chose 2-layer snow estimates to match the observations. The reason for the MCMC choice and the possible influence from the prior are discussed. Analysis of Response of Microwave TB to SWE To better understand the reason for the outliers and the underestimation for the CLPX snowpits, Figure 29 shows the measured and the synthetic TB versus SWE for all the 48 snowpits at different frequencies. It can be found that, both the measured TB and the synthetic TB at 36.5 GHz show a decreasing trend as the SWE increases from 40 mm to about 120 mm. This trend becomes unclear for SWE greater than 120 mm due to the influence of other snow parameters. This explains why the underestimation for CLPX snowpits can happen because their SWE is inside this insensitive region. More information can be found in simulation test in Figure 30. The measured multiple-layer snow density, temperature and pec were used for simulation for each snowpit, but the total snow depth was scaled to match the different SWE in X-axis. Each line corresponds to a snowpit, and the relative thickness of layers was not changed. In general, it shows that, the simulated TB at 36.5 GHz decreases with the increasing SWE from 20 mm up to about 160 mm. After 160 mm, the simulated TB begins to increase with increasing SWE, 98 for some snowpits. This trend is also found in the field measurements for tundra snow in Canada (Derksen et al., 2010) and alpine snow in Switzerland (Matzler, 1982). The inflection point from decreasing to increasing trend is determined by snow properties, especially, the snow micro-structure parameter. This complex TB response indicates that it can be hard to do SWE retrieval for deep snow, unless better information of other snow properties is available. Another interesting thing is the sensitivity of 89-GHz TB to SWE, as can be found both in measurements (Figure 29(a)) and simulations (Figure 29(b) and Figure 30). The reason for the sensitivity will be discussed later in the next subsection. For some snowpits, the 89-GHz TB is sensitive to SWE when the 36.5-GHz TB is not. It indicates the inclusion of this frequency may help improve the SWE retrieval if it is used properly. For the lower frequencies, the results based on the both measured and synthetic TB in Figure 29 show that TB at 6.925 and 10.65 GHz is insensitive to SWE. This is because, according to RT theory, the dimension of snow particles is much smaller than the long wavelength at such frequencies, and the volume scattering effect is weak; these channels may aid in estimation of soil properties. In Figure 29(a), there are two snowpits with the measured TB at 18.7 GHz lower than 240 K, which is not found for snowpits with similar SWE or in their corresponding synthetically-generated TB in Figure 29(b). These are the two Churchill snowpits where the strongly overestimated BASE-PM SWE was produced. Figure 30 demonstrates why such overestimation is possible, because in general, with long wavelength, TB at 18.7 GHz is only slightly sensitive to SWE for shallow snow; this sensitivity becomes more apparent only for very deep snow. According to simulation test, 99 the very low TB at 18.7 GHz but not 10.65 GHz may be related to the coherent effects of ice lenses with some specific thicknesses. It indicates, although 18.7-GHz remains sensitive for deep snow, other factors also need to be considered. For these two snowpits, the retrieved SWE can be improved by excluding this frequency. (b) Simulated TB 280 280 260 260 240 240 220 220 200 200 TB (K) TB (K) (a) Measured TB 180 160 160 6.925GHz 10.65GHz 18.7GHz 36.5GHz 89GHz 140 120 100 0 180 50 100 150 SWE (mm) 200 6.925GHz 10.65GHz 18.7GHz 36.5GHz 89GHz 140 120 100 0 250 50 100 150 SWE (mm) 200 250 Figure 29. TB at 6.925 to 89 GHz versus SWE for all the 48 snowpits, where in (a) it shows the measured TB, in (b) it shows the synthetic TB. 100 (b) (c) 280 260 260 260 240 240 240 220 220 220 180 200 B 200 T (K) 280 TB (K) B T (K) (a) 280 180 200 180 160 160 160 140 140 140 120 120 120 100 0 100 200 300 SWE (mm) 400 100 0 100 200 300 SWE (mm) 400 100 0 10.65GHz 18.7GHz 36.5GHz 89GHz 100 200 300 SWE (mm) 400 Figure 30. Simulated TB at 10.65 to 89 GHz using the measured snow density, temperature and pec for each snowpit but the snow depth was scaled to match the SWE from 20 mm to 400 mm in X-axis. Each line corresponds to the TB-SWE relationship for one snowpit at one frequency: (a) for the snowpits with SWE between 40 and 50 mm; (b) for the snowpits with SWE between 50 and 170 mm; (c) for the snowpits with SWE larger than 170 mm (i.e. the CLPX snowpits). The Influence of The Number of Snow Layers The issues related to the number of snow layers (nlyr) can be summarized in the following three questions: first, how many snow layers does the snowpit truly have? Second, how many snow layers are required to reproduce the measured TB signal? Finally, as mentioned in Section 4.4, will the prior for the number of snow layers influence the retrieved SWE? The natural snow cover is vertically-inhomogeneous and could have many layers. For the Sodankylä snowpits, usually 5 to 10 layers were recorded. The Churchill snowpits measured by the Environment Canada group had 4 to 7 layers. The most high-frequent number of layers is 2 for the entire CLPX study region (Durand and Liu, 2012), and 5 for the Fraser MSA. 101 In this study, due to the computational cost, BASE-PM runs with a maximum of 6 layers. In most cases, the algorithm chose 2 layers. Using the generic prior, more layers were not chosen, because it was found that the observed TB at all frequencies utilized can be reproduced well by a 2-layer profile. On the other hand, it was found at least two layers were required for some snowpits to match the TB observations, especially when the measured TB at 89 GHz is higher than 36.5 GHz. Figure 31 shows the SWE estimates when BASE-PM was forced to choose 1 layer. The rms error of SWE is 79.2 mm for all snowpits, and 73.2 mm excluding the two Churchill snowpits. For measured SWE between 90 to 250 mm, the estimated SWE for most snowpits is between 160 and 190 mm. Moreover, the uncertainty of the SWE estimates is dramatically underestimated. Presumably, this is because constraining the snowpack to a single layer essentially presupposes a simpler radiative transfer dynamic that does not hold for these snowpits. 102 (a) (b) MCMC SWE uncertainty (mm) SWEMCMC(mm) 400 300 sd−nov sd−dec sd−jan sd−feb sd−mar ch−jan ch−feb ch−mar ch−apr clpx−feb 200 100 0 0 100 200 300 SWEtrue (mm) 400 300 250 200 150 100 50 0 0 100 200 |SWE error| (mm) 300 Figure 31. The SWE using 1-layer snow parameter estimation when the same priors as BASE-PM were used: (a) the MCMC SWE versus the true SWE; (b) the estimated SWE uncertainty versus the absolute SWE error. To understand the reasons, Figure 32(a) shows the simulated TB versus SWE for all snowpits like in Figure 30, but the simulations here used mass-weighted (1-layer) snow properties in the first row, and 2-layer snow properties (relayered from the original snowpit measurements) in the second row. It shows when snow has only 1 layer, the simulated TB at 89 GHz will not change when SWE is larger than 40 mm. However, it conflicts with the true SWE-TB relationship according to measurements in Figure 29 and simulations using originally-measured multiple-layer snow profiles in Figure 30. The increasing trend of 89-GHz TB with SWE can be reproduced using 2-layer snow. The reason for the sensitivity of 89-GHz to SWE is related to stratigraphy. 103 (a) 260 240 240 240 220 220 220 200 200 200 180 B B T (K) 260 180 180 160 160 160 140 140 140 120 120 120 100 200 300 SWE (mm) (d) 280 100 200 300 SWE (mm) 100 0 400 (e) 280 280 260 260 240 240 240 220 220 220 200 200 200 180 B B 180 T (K) 260 T (K) B 100 0 400 160 160 140 140 140 120 120 120 100 200 300 SWE (mm) 400 100 0 100 200 300 SWE (mm) 400 10.65GHz 18.7GHz 36.5GHz 89GHz 100 200 300 SWE (mm) (f) 400 180 160 100 0 (c) 280 260 100 0 T (K) (b) 280 T (K) TB (K) 280 100 0 10.65GHz 18.7GHz 36.5GHz 89GHz 100 200 300 SWE (mm) 400 Figure 32. Simulated TB at 10.65 to 89 GHz using the 1-layer (first row) or 2-layer (second row) snow density, temperature and pec for each snowpit, where the snow depth was scaled to match the SWE from 20 mm to 400 mm in X-axis. Each line corresponds to the TB-SWE relationship for one snowpit at one frequency. From left to right: the first column is for the snowpits with SWE between 40 and 50 mm; the second column is for the snowpits with SWE between 50 and 170 mm; the third column is for the snowpits with SWE larger than 170 mm (i.e. the CLPX snowpits). 104 (a) (b) 0 -5 -15 -25 -25 -25 -35 -35 -35 -55 z (cm) -15 -45 -45 -55 -45 -55 -65 -65 -65 -75 -75 -75 -85 140 180 220 T B (K) 260 300 -85 140 180 220 T B (K) (c) 0 -5 -15 z (cm) z (cm) 0 -5 260 300 -85 140 180 220 T B (K) 260 300 Figure 33. Simulated upward propagated TB at 89 GHz, vertical polarization for the example snowpit. From left to right, 1-layer, 2-layer and the originally-measured 6-layer snow profiles were utilized. The originally-measured snow depth was used for the black line; for blue dash line and red line, the snow depth was scaled by 2/3 and 1/3, respectively, with relative thickness of layers unchanged. Figure 33 shows the upward-propagated TB at 89 GHz inside snowpack using the method in Pan et al. (2016) when the 1-layer, 2-layer and the originally-measured 6-layer snow profiles are used for the example snowpit. This snowpit has two surface layers with a pec of 0.09 mm and a density smaller than 200 kg/m3. Figure 33 shows that the decrease of upward-propagated TB is slower in the surface layers than the bottom layers. The changing rate of upward-propagated TB is more influenced by pec than density. When there is only 1 layer, no matter if the snow depth is 84 cm, or 2/3 or 1/3 of the original depth, TB decreases to the same as about 142 K at the snow surface. However, when there are at least two layers, the snow profile with thicker snow depth has thicker surface layer with small volume scattering. The scattering coefficient is smaller than the absorption 105 coefficient. In this case, according to radiative transfer theory, the snow-emitted TB in this layer increases the TB observed at the snow surface. The simulation here assumed unchanged relative thickness of layers. However, if the density and thickness of the surface layer are kept unchanged and only the thickness of the bottom layer is scaled to match the 2/3 or 1/3 SWE, MELMS-IBA gives a similar TB (200.8 K) as the original snow profile. Therefore, it indicates that the increase in 89-GHz TB with increasing SWE is related to the thickness of the surface layer with smaller pec. In natural conditions, the snow compaction and grain growth are gradual and continuous processes throughout the entire snowpack. It is difficult to keep the thickness of surface layer unchanged, but increasing only the SWE amount in the bottom layers. Therefore, sensitivity of 89-GHz TB is partly related to the snow emission theory, but also influenced by the natural properties of snow. Additionally, related to the effects of layered snow, note that TB at 36.5 GHz will increase with increasing SWE when SWE is larger than a certain value using the measured snow profiles, but this feature can also not be captured by a 1-layer model. Therefore, considering at least two snow layers is important for passive microwave SWE estimation, especially for deep snow case. To test the influence of the prior for number of snow layers (nlyr), SWE retrieval using different λ for Poisson distribution or a uniform distribution were tested. Results show that the algorithm continued to choose 2 layers as long as generic stratigraphy prior was used. When local stratigraphy prior was used, more layers can be chosen for some 106 snowpits, but the resulted change in SWE for these snowpits was only between 2 to 8 mm. Therefore, in general, the prior for nlyr is not important. 107 Chapter 4: SWE Retrieval Using Active Microwave Measurements In this chapter, the BASE algorithm using scatterometer observations and generic priors is tested for the Sodankyla SnowScat and snowpit measurements (Lemmetyinen et al., 2016). IOP1 is taken as an example due to the computational cost. Section 1 introduces the observations utilized and the estimated parameters. Section 2 and 3 show the retrieval results using the real SnowScat measurements and synthetic-generated backscattering coefficients, respectively. Section 4 discusses the algorithm performance. Section 4.1: Active Microwave Retrieval Algorithm Setup Section 2.3 shows MEMLS-IBA can simulate multi-frequency TB with unbiased sensitivity to snow depth and snow grain size. Section 2.4 shows MEMLS3&a with IBA can simulate the SnowScat backscattering coefficient measurements with acceptable accuracy when the snowpacks are not influenced by partly frozen soil. Therefore, here MEMLS3&a with IBA was chosen as the observation model for the active microwave SWE retrieval. The estimated parameters include the multi-layer snow thickness, temperature, density and exponential correlation length (pec), and soil moisture, temperature and roughness, similar to the passive microwave SWE retrieval algorithm. The polarization splitting parameter qs was also estimated, with a prior as 0.1±0.02; m is fixed at 0.1. The measured 108 backscattering coefficient was assumed to have a normal-distributed error with 0.1-dB standard deviation. Section 4.2: Retrieval Results Using SnowScat Measurements The backscattering coefficients at 10.2, 13.3 and 16.7 GHz, VV-polarization observed in 50° incidence angle at 2 AM every other day were extracted and input into the MCMC algorithm. It should be noted that the VH-pol. backscattering coefficient measurements were not used, because the retrieval test using 3-frequency, VV- and VH-pol. SnowScat measurements did not give reasonable snow estimates. The MCMC algorithm favored very large pec (above 0.6 mm) and very small snow depth when cross-pol observations were used. The reason is most likely due to the errors in MEMLS3&a or the errors in the SnowScat measurements; in other words, the model can not reproduce all the SnowScat measurements with a same level of unbiased errors using reasonable values of snow parameters. Therefore, the VH-pol. measurements are not used when the algorithm uses 3-frequency measurements. Figure 34 shows an example of the Markov chains retrieved from the SnowScat measurements on February 4th, 2010. The chains of the layer thickness, pec and snow density of the bottom and the surface layers are shown on the left; the histograms of the chains and the priors are shown on the right. Figure 35 shows that the backscatter coefficient simulated from the MCMC estimates matches the SnowScat observations. Note that the chains have 300,000 iterations, which are 10 times of the iterations utilized for passive microwave SWE retrieval. The much larger chain length is because of the low 109 convergence speed. Indeed, the first 30,000 iterations (i.e., the first 10% of chains) sampled only the relatively large pec and low layer thickness values of the posterior MCMC histograms. Chain convergence is required to ensure that the posterior actually satisfies Bayes law. Secondly, it can be seen from Figure 34 that, the posterior histograms of layer thickness, pec and density for the surface layer (Lyr2) are very close to their prior histograms. It means only the bottom-layer snow properties were constrained by the backscattering coefficient measurements. 110 Figure 34. Markov chains of layer thickness (dz), exponential correlation length (pec) and snow density for the bottom (Lyr1) and the surface (Lyr2) snow layers (left), and their posterior histogram compared to the prior histogram (right). 111 Figure 35. Backscattering coefficient calculated from the snow parameters in Markov chains in Figure 34. Figure 36 shows the sensitivity of the backscattering coefficient to thickness, pec and density at the bottom and the surface snow layers for the example in Figure 34. The simulations here are based on the 2-layer MCMC snow estimates. In each figure, the cross indicates the original MCMC estimate for the parameter of interest; this parameter varies with other parameters unchanged to produce the corresponding changes in backscattering coefficients (dotted lines). It can be seen from Figure 36 that, in general, the backscattering coefficient is not very sensitive to the properties of the surface snow layer, especially the surface-layer density and thickness. It explains their similarities of posterior histogram to the prior histogram in Figure 34. Also, it can be found that, the backscattering coefficient is both sensitive to the bottom-layer thickness and pec. However, it seems the sensitivity is stronger for small bottom-layer thickness, but for 112 larger bottom-layer pec. The backscattering coefficient is also not very sensitive to the snow density. The sensitivity is slightly higher for the bottom-layer because the bottomlayer density determines the refraction at the snow-soil boundary. Snow density is involved in the snow scattering coefficient calculation in IBA, and its influence can be slight larger for the bottom layer with a large pec. 113 Bottom-Layer: sensitivity to Dz -5 -7 -9 -9 -11 -11 -13 -13 VV (dB) 0 VV (dB) 0 -7 -15 -15 -17 -17 -19 -19 -21 0 50 100 150 200 layer thickness (dz) (cm) -25 -9 -13 VV (dB) 0 -13 -15 -17 -19 -19 -21 -21 -23 -23 0.2 0.3 0.4 p ec (mm) -25 -9 -13 VV (dB) 0 -13 -15 -17 -19 -21 -23 -23 600 3 snow density (kg/m ) 0.2 0.3 0.4 p ec (mm) 800 0.5 Surface-Layer: sensitivity to density 10.2VV 13.3VV 16.7VV -17 -21 400 0.1 -15 -19 200 0 -9 -11 0 10.2VV 13.3VV 16.7VV -7 -11 -25 200 Surface-Layer: sensitivity to Pec -5 10.2VV 13.3VV 16.7VV -7 VV (dB) 0 0.5 Bottom-Layer: sensitivity to density -5 150 -15 -17 0.1 100 -9 -11 0 50 layer thickness (dz) (cm) -7 -11 -25 0 -5 10.2VV 13.3VV 16.7VV -7 10.2VV 13.3VV 16.7VV -23 Bottom-Layer: sensitivity to Pec -5 VV (dB) 0 -21 10.2VV 13.3VV 16.7VV -23 -25 Surface-Layer: sensitivity to Dz -5 -25 0 200 400 600 snow density (kg/m 3 ) 800 Figure 36. Sensitivity of backscattering coefficient to the bottom- and surface-layer thickness, pec and density for the example retrieval results in Figure 34. 114 Figure 37 shows the MCMC-estimated snow depth, density and SWE for all the SnowScat measurements utilized for retrieval in IOP1, based on 2-layer snow estimates. The retrieved snow depth is compared with the snow depth measured by the acoustic Campbell SR50 sensor, and the SWE is compared with the Gamma Water Instrument (GWI) measurements. The depth-averaged density is compared with the depth-averaged density from the closest snowpit measurements in time. Basically, the retrieval results can be classified in three parts. Before January 27, the snow depth is estimated with good accuracy, but the retrieved snow density biased high for about 50%, and it resulted in an overestimated SWE for about 70-mm. After February 15, the retrieval algorithm gives smaller snow depth compared to the measurements, and results in an underestimated SWE. Figure 38 takes point 1 in Figure 37 as an example to show the reasons for the SWE overestimation. It can be found that, the retrieval algorithm estimated a very large bottom-layer snow density to match the observations at three frequencies. Figure 39 takes point 29 as an example to show the reason for the SWE underestimation. It can be found that, the retrieval algorithm accurately estimated bottom-layer thickness and pec; however, it underestimated the surface-layer thickness. For all the points after January 27, it was found that the surface-layer thickness is very close to the prior thickness. The underestimated snow depth is mainly because of the underestimated surface-layer thickness, which is insensitive to the backscattering coefficients at the frequencies utilized here. 115 IOP1 - Retrival Using Real Observation 90 80 measured mcmc prior 70 SD (cm) 12 60 3 50 1 4 2 17 14 56 7 8 22 18 16 27 29 24 21 15 13 19 20 33 32 34 23 11 9 31 28 30 26 10 40 25 30 Jan Feb 350 6 1 Snow density (kg/m 3 ) 2 3 4 Mar Apr 7 8 measured mcmc prior 5 9 300 11 12 17 18 13 250 19 10 15 14 16 28 20 21 22 23 25 Feb 32 34 33 26 24 200 Jan 27 29 31 30 Mar Apr 200 180 3 6 7 2 160 4 1 5 12 8 17 SWE (mm) 9 13 140 11 120 31 18 15 16 14 28 30 19 26 22 20 29 32 34 33 27 24 10 21 100 80 60 Jan 23 measured mcmc prior 25 Feb Mar Apr Figure 37. MCMC-retrieved snow depth, snow density (bulk values) and SWE (from top to bottom) using 3-frequency SnowScat observations at VV pol. based on 2-layer snow estimates 116 50 pit1: density MCMC True 40 35 35 30 30 30 z (cm) 40 35 25 25 20 20 20 15 15 15 10 10 10 5 5 5 0 100 0 200 300 400 0 3 0.2 MCMC True 45 40 25 pit1: temperature 50 MCMC True 45 z (cm) z (cm) 45 pit1: micro-structure 50 0 -15 0.4 -10 density (kg/m ) -5 0 T (o C) p ec (mm) Figure 38. Retrieved snow profile for point 1 in Figure 34. 80 pit29: density MCMC True 80 MCMC True 70 60 50 50 50 z (cm) 60 40 40 30 30 20 20 20 10 10 10 0 200 250 3 density (kg/m ) 300 0 0.2 MCMC True 40 30 0 150 pit29: temperature 70 60 z (cm) z (cm) 70 pit29: micro-structure 80 0.4 0 -15 -10 -5 0 T (o C) p ec (mm) Figure 39. Retrieved snow profile for point 29 in Figure 34. Figure 40 shows the backscattering coefficients calculated from the MCMC estimates and their errors. It can be found that the errors are smaller than the MEMLS3&a simulation errors using snowpit measurements in Figure 14, because the retrieval algorithm tried its best to tune all the estimated snow and soil parameters to match the observations. On the other hand, there is still an error of about 0.1 to 0.6 dB, and the bias is systematic at different frequencies. It indicates the deficiency of the observation model 117 to explain the SnowScat measurements. In Figure 14, the MEMLS3&a simulation underestimated the backscattering coefficient at 10.2 and 16.7 GHz in January and it underestimated the backscattering coefficient at 16.7 GHz in March. More importantly, the features of the biases in Figure 14 are different with Figure 40. Therefore, the corresponding snow properties in the MCMC retrieval and in the snowpit measurements can be different. -8 MCMC-Chain Simulations from MCMC snow/soil estimates -10 -12 vv 0 -14 -16 -18 10.2VV-Obsr 10.2VV-MCMC 13.3VV-Obsr 13.3VV-MCMC 16.7VV-Obsr 16.7VV-MCMC -20 -22 Jan 0.8 Feb Mar Apr Error of MCMC-Chain Simulations 0.6 0.4 10.2VV-error 13.3VV-error 16.7VV-error vv 0 0.2 0 -0.2 -0.4 -0.6 Jan Feb Mar Apr Figure 40. Backscattering coefficient simulated from the MCMC estimates (top) and their errors compared to the SnowScat measurements (bottom). 118 90 80 measured mcmc prior SD (cm) 70 60 50 40 30 Jan Feb Mar 380 measured mcmc prior 360 Snow density (kg/m 3 ) Apr 340 320 300 280 260 240 220 200 Jan Feb Mar Apr IOP1 200 180 SWE (mm) 160 140 120 100 measured mcmc prior 80 60 Jan Feb Mar Apr Figure 41. MCMC-retrieved snow depth, snow density (bulk values) and SWE (from top to bottom) using 3-frequency SnowScat observations at VV pol. based on 1-layer snow estimates. Figure 41 shows the estimated snow depth, snow density and SWE when the retrieval algorithm chooses 1-layer snow estimates. It can be found that the retrieval results using 1-layer snow estimates are very close to the results using 2-layer snow estimate with only 119 slight differences. It also overestimates the density for shallow snow in January and underestimates the depth for deep snow in March. In summary, the retrieved SWE using the actual SnowScat measurements has low correlation with the measured SWE. Using 2-layer estimates, the mean bias of the MCMC-estimated snow depth is -9.8 cm, and the rms error is 13.5 cm. The mean bias of SWE is 2.5 mm, and the rms error is 44.2 mm. The algorithm overestimates SWE when the snow density is overestimated for shallow snow, and underestimates SWE when the surface-layer thickness insensitive to the backscattering coefficient is underestimated for deep snow. Section 4.3: Retrieval Results Using Synthetic Backscattering Coefficient In the synthetic test, MEMLS3&a uses the measured snow depth from the acoustic sensor at the same time of the backscattering measurements to generate the synthetic backscattering coefficient values. The snowpit measurements were relayered into 2-layer values to provide the snow stratigraphy information as the snow density, temperature and pec. The snow depth of the snowpit was scaled to match the snow depth from acoustic sensor, with the relative thicknesses of snow layers unchanged. The soil temperature and moisture are from the 2-cm Delta-T ML2x and Pentronic PT100 sensor measurements (Lemmetyinen et al., 2016). The soil roughness was fixed as 0.1 cm. The “true” SWE is calculated from the 2-layer snow thickness and density. Figure 42 shows the retrieved SWE using synthetically-generated backscattering measurements when the algorithm estimated 2-layer snow properties. The results include 120 34 points. The estimated snow density is about 217 kg/m3 for all points (not shown here), which reproduced the density priors. Therefore, the error in SWE is mainly from the errors in the layer thickness estimates. The mean bias of the retrieved SWE is -30.8 mm and the rms error is 42.2 mm; it improves over the SWE prior, which has a mean bias of 53.1 mm and a rms error of 58.1 mm. The mean bias of the retrieved snow depth is -9.4 cm, and the rms error is 12.5 cm; it also improves over the snow depth prior, which has a mean bias of -23.2 cm and a rms error of 23.6 cm. IOP1 - Synthetic Test - 2-layer 250 SWE (mm) 200 mcmc prior measured 26 17 150 13 6 100 12 3 45 7 8 9 50 Jan 16 20 21 24 27 30 32 33 25 15 Feb 34 29 31 18 12 14 11 10 28 19 2223 Mar Apr Figure 42. MCMC-retrieved SWE using synthetic-generated backscattering coefficients based on 2-layer snow estimates. Figure 43 compares the 2-layer MCMC-estimated and true snow thickness and pec. The true values are exactly the values used in MEMLS3&a to generate the backscattering coefficient observations. It can be seen that, the estimated bottom-layer pec has high correlation with the true pec. However, the surface-layer thickness closely matches the prior thickness. It agrees with the small sensitivity of this parameter shown in Figure 36. 121 The MCMC-estimate bottom- and surface-layer pec are higher than the true values. It is mostly likely compensated by the lower bottom-layer thickness estimated by MCMC. 0.4 Bottom Layer 0.35 Pec (mm) 0.2 0.05 Jan 0.8 Layer thickness (m) 0.7 0.2 0.15 Obsr MCMC Prior Feb 0.1 Mar Apr 0.05 Jan Bottom Layer 0.7 Obsr MCMC Prior 0.6 0.6 0.5 0.4 0.3 0.2 0.1 Jan 0.25 Layer thickness (m) Pec (mm) 0.25 0.1 Obsr MCMC Prior 0.3 0.3 0.15 Surface Layer 0.35 Feb Mar Apr Surface Layer Obsr MCMC Prior 0.5 0.4 0.3 0.2 0.1 Feb Mar Apr 0 Jan Feb Mar Apr Figure 43. Comparison of retrieved and true bottom- and surface-layer pec and thickness based on 2-layer snow estimates. Figure 44 shows the retrieved SWE using synthetically-generated backscattering coefficients based on 1-layer snow estimates. The results in Figure 44 are better than the results using 2-layer snow estimates in Figure 42. The mean bias and the rms error for SWE are -7.5 mm and 34.1 mm, respectively; and the mean bias and the rms error for snow depth are -0.9 cm and 11.1 cm. Statistically, the SWE rms error using 1-layer snow 122 estimates is close to the passive microwave SWE retrieval accuracy in Chapter 3. However, the active microwave retrieval experiment here utilized many more iterations, and it only worked for the synthetically generated microwave measurements. IOP1- Synthetic Test - 1-layer 250 SWE (mm) 200 mcmc prior measured 150 100 50 Jan Feb Mar Apr Figure 44. MCMC-retrieved SWE using synthetic-generated backscattering coefficients based on 1-layer snow estimates. In Figure 42, the retrieved SWE is significantly underestimated for point 9-11 and 19-25 using 2-layer snow estimates. In Figure 45 and 46, point 9 and 22 are taken as examples to explore the reasons for the errors. From Figure 45, it can be found that the error for point 9 using 2-layer snow estimates is because the input snow profiles have basically homogenous snow properties except a very thin surface layer. Therefore, the retrieved SWE using 1-layer snow estimates is much closer to the SWE observations. The situation is similar the point 22 in Figure 46. However, the estimated SWE is also not very good using 1-layer snow estimates. This may be because the pec for this snowpit is too small. The total magnitude of volume scattering required to match the observations is so small 123 that a small error in estimated pec can result in a large error in snow depth estimate, and then lead to the error in SWE. pt9: density pt9: micro-structure 60 50 40 40 40 30 10 0 100 z (cm) 50 20 30 20 10 MCMC True 200 250 60 30 10 MCMC True 0 150 0 0.2 0 -10 0.4 pt9: micro-structure 60 60 50 40 40 40 20 10 0 100 z (cm) 50 z (cm) 50 30 30 20 200 density (kg/m 3 ) 250 -6 -4 pt9: temperature MCMC True 30 10 MCMC True 0 150 -8 20 10 MCMC True MCMC True T (o C) p ec (mm) pt9: density pt9: temperature 20 density (kg/m 3 ) z (cm) 60 50 z (cm) z (cm) 60 0 0.2 0.4 0 -15 -10 -5 0 T (o C) p ec (mm) Figure 45. Retrieved 1- (top) and 2-layer (bottom) snow profiles for point 9 based on the synthetic-generated backscattering coefficients in Figure 42 and 44. 124 pt22: density 70 pt22: micro-structure 50 50 50 40 40 40 30 10 0 100 MCMC True 150 200 30 20 10 10 0 0.05 250 0 -30 0.15 pt22: density 70 pt22: micro-structure 70 60 50 50 50 40 40 40 20 10 0 100 MCMC True 150 200 250 30 20 10 10 density (kg/m 3 ) 0.1 0 0.15 pt22: temperature MCMC True 30 20 0 0.05 -10 60 MCMC True z (cm) z (cm) 60 30 -20 T (o C) p ec (mm) density (kg/m ) 70 0.1 MCMC True 30 20 3 pt22: temperature 60 z (cm) 60 20 z (cm) 70 MCMC True 60 z (cm) z (cm) 70 0 -30 -20 -10 0 T (o C) p ec (mm) Figure 46. Retrieved 1- (top) and 2-layer (bottom) snow profiles for point 22 in Figure 42 and 44. In Figure 47, the snow profiles for point 30 explain the reasons for the slightly underestimated SWE using 2-layer snow estimates for the points 1-4 and 29-33. It shows the SWE error is mainly from the bias in the density estimates. The estimated snow density (215.5 kg/m3) at the bottom layer reproduced the prior (217 kg/m3), and is smaller than the true snow density (239.3 kg/m3). The underestimated bottom-layer density, together with the underestimated surface layer thickness give a 30-mm biased-low SWE compared to the true SWE. 125 pt30: density pt30: micro-structure 80 70 60 60 60 50 50 50 40 0 150 40 30 20 10 z (cm) 70 30 20 250 0 pt30: density 0.2 0 -10 0.3 80 MCMC True 70 50 z (cm) 60 50 z (cm) 60 50 40 40 30 30 20 20 20 10 10 10 0 200 250 0 0.2 -4 pt30: temperature 40 30 density (kg/m 3 ) -6 70 60 0 150 -8 T (o C) pt30: micro-structure 80 MCMC True 70 0.1 p ec (mm) density (kg/m ) 80 40 10 0 200 MCMC True 20 MCMC True 10 MCMC True pt30: temperature 30 3 z (cm) 80 70 z (cm) z (cm) 80 0.4 0 -15 MCMC True -10 -5 T (o C) p ec (mm) Figure 47. Retrieved 1- (top) and 2-layer (bottom) snow profiles for point 30 in Figure 42 and 44. Section 4.4: Discussion Two major problems are found for SWE retrieval using active backscattering measurements: the overestimated SWE in January due to the overestimated bottom-layer snow density, and the underestimated SWE in March due to the insensitivity of surface layer thickness to the backscattering coefficients. For the points in January, the overestimated density and SWE are not found in the synthetic test. In the synthetic test, the posterior snow density is determined by the prior. It implies that the algorithm will not tune snow density if other snow and soil parameters can reproduce the backscattering coefficient observations. The MEMLS3&a simulation 126 results using the snowpit measurements (in Figure 14) show that MEMLS3&a underestimates the backscattering coefficient at 10.2 and 16.7 GHz with an unbiased 13.3-GHz simulation. The sensitivity test for point 1 (see Figure 48) shows an increased bottom-layer thickness or pec increases the backscattering coefficients at three frequencies and increases the frequency differences. However, changing the bottom-layer density from about 230 kg/m3 (measurements) to 380 kg/m3 (MCMC estimates) decreases the backscattering coefficient at 10.2 GHz when it increases the frequency differences. Also, while the influences of the layer thickness and pec to backscattering coefficient are stronger at 16.7 GHz, the influence of snow density is stronger at 10.2 GHz. These three parameters mutually influence each other to reproduce the 3-frequency observations. For point 1, the influence of the bottom-layer density cannot be compensated by any other parameters, including the soil parameters. The fact that the overestimation is not found in the synthetic test indicates the overestimated bottom-layer density is related to the special pattern of errors in the observation model compared to the actual SnowScat measurements. Such model error may be in the snow RT model or the soil reflectivity model, because the errors here is found for shallow snow (<100 mm). Therefore, it suggests an improved observation model to improve the SWE retrieval accuracy. 127 Bottom-Layer: sensitivity to Dz 5 10.2VV 13.3VV 16.7VV -5 10.2VV 13.3VV 16.7VV 0 VV (dB) 0 VV (dB) 0 0 Surface-Layer: sensitivity to Dz 5 -10 -5 -10 -15 -15 -20 -20 -25 -25 0 50 100 150 200 0 layer thickness (cm) 5 -10 200 10.2VV 13.3VV 16.7VV -5 -10 -15 -15 -20 -20 -25 -25 0 0.1 0.2 0.3 0.4 0.5 0 0.1 layer thickness (dz) Bottom-Layer: sensitivity to density 0.4 0.5 Surface-Layer: sensitivity to density 10.2VV 13.3VV 16.7VV 0 VV (dB) 0 -5 0.3 5 10.2VV 13.3VV 16.7VV 0 0.2 p ec (mm) 5 VV (dB) 0 150 Surface-Layer: sensitivity to Pec 0 VV (dB) 0 VV (dB) 0 5 10.2VV 13.3VV 16.7VV -5 100 layer thickness (cm) Bottom-Layer: sensitivity to Pec 0 50 -10 -5 -10 -15 -15 -20 -20 -25 -25 0 200 400 600 snow density (kg/m 3 ) 800 0 200 400 600 800 snow density (kg/m 3 ) Figure 48. Sensitivity of backscattering coefficient to the bottom- and surface-layer thickness, pec and density for point 1 in Figure 37. The underestimated SWE in March using the actual SnowScat measurements, as analyzed above, may be influenced by the insensitivity of the backscattering coefficient to surface-layer thickness. However, the fact that a strong underestimation is found also for the real test based on 1-layer snow estimates indicates that it is likely also influenced by the observation model error. In Figure 14, MELMS3&a gives an underestimated backscattering coefficient at 16.7 GHz but the simulations are unbiased at 10.2 and 13.3 128 GHz. The MCMC algorithm tuned pec larger and the snow depth smaller to remove these frequency-dependent errors, because the “observation error” given to MCMC is the same at all frequencies. It then produces the errors in the MCMC retrieval results. On the other hand, it is true that the surface-layer thickness is not very sensitive to the backscattering coefficients from 10.2 to 16.7 GHz. Large errors can be produced for snowpacks with a very thick fresh snow layer and a thin layer of melting crusts from the last snowfall event. As the surface layer thickness is strongly influenced by the prior, it may be possible to reduce SWE error by providing a good prior for the surface layer thickness, for example, by using some prior for the relative thicknesses of the surface and the bottom layers. Also, the passive microwave SWE retrieval study mentioned the sensitivity of 90-GHz TB to the surface layer thickness. Therefore, another option is to use a synergy of radar with the TB measurements at higher frequencies. In both the real and the synthetic tests, the estimated SWE based on 1-layer snow estimates is very close to, or sometimes even improves over the SWE based on 2-layer snow estimates. It indicates that maybe only a single layer is required for the active microwave SWE retrieval instead of multiple layers. This may be because the penetration depth is large for frequencies from 10.2 to 16.7 GHz. When the underlying soil is homogenous, it will not produce strong differences in the total magnitude of snow volume scattering using one or multiple snow layers; this assumption is supported by the forward model simulations in Du et al. (2010b). It should be noted that, future studies using a different snow RT model or a multi-layer frozen soil model may give a different conclusion. 129 Chapter 5: Conclusions SWE estimation is of critical importance for hydrological applications. This research developed a new SWE retrieval algorithm based on a Bayesian-based approach to iteratively find the posterior probability of snow and soil parameters constrained by the multi-frequency passive or active microwave measurements. To explore appropriate snow radiative transfer (RT) models that can be used as forward snow microwave observation models, two commonly-used models, the Microwave Emission Model of Layered Snowpacks (MEMLS) and the Helsinki University of Technology (HUT) snow emission model, were compared and evaluated in theory and in practice. Three major theoretical differences were identified: the 1-flux versus 2-flux RT theories, the different volume scattering coefficients, and the trapped radiation that is not considered in the HUT model. It was found that the HUT model tends to underestimate TB at the snow surface for thick snow due to the simplification of RT equation resolution. MEMLS with the Improved Born approximation (IBA) performed best among the RT equation resolution and the volume scattering models. It gives an rms error of about 10 K for TB simulation, and an rms error of about 1 dB for backscattering coefficient simulation. MELMS3&a (Microwave Emission Model of Layered Snowpacks 3&active) 130 underestimated the backscattering coefficient for snowpacks with partly-frozen underlying soil conditions. When MEMLS-IBA is used as the observation model in the Markov Chain Monte Carlo (MCMC) retrieval system, it achieved an accuracy of about 30.8 mm for shallow snow (excluding deeper snow and two outliers). The algorithm is applicable requiring only a set of globally-available prior information. The estimation of 2-layer snow properties is suggested to match all the TB measurements from 10.65 to 90 GHz. TB at 90 GHz is sensitive to the thickness of surface snow layers with small grain size and volume scattering. When the retrieval algorithm is adapted for the active microwave SWE retrieval, however, the retrieval was successful only in the synthetic test. A similar SWE retrieval accuracy is found utilizing the synthetically-generated backscattering coefficients and based on 1layer snow estimates. Also, it was found the backscattering coefficient from 10.2 to 16.7 GHz is not sensitive to the thickness of the surface snow layers. Such layers are “transparent” for X- and Ku-band radar radiation. The retrieval result suggests an improvement in the active snow RT model or the soil RT model. The small sensitivity of surface-layer snow properties to the active microwave measurements suggests a better prior or a synergy of active and passive microwave measurements at high frequencies. Continued improvement in the understanding of microwave signals and microwave simulation models is important to support the satellite SWE retrieval algorithms, especially in regions with a sparse conventional observing network, such as the boreal forest and (sub) Arctic tundra. The benefit of the retrieval algorithm developed in this 131 research is that it can be used for different snow types in different regions due to the physically-based observation model it utilized. It does not require ground measurements or land surface model simulations, although they may be used to provide better prior information; users can employ application-specific configurations. It can incorporate the measurements from different sources, as long as the observation models have adequate fidelity. Therefore, the framework can be expanded to types of remotely-sensed measurements. The algorithm can give an uncertainty estimate for its retrieved SWE, which can help recognize outliers without the ground validation datasets. The uncertainty estimates will be helpful in a data assimilation framework when the retrieved SWE is merged with the estimation from other techniques. 132 Bibliography Andreadis, K. M., Lettenmaier, D. P. Implications of representing snowpack stratigraphy for the assimilation of passive microwave satellite observations. Journal of Hydrometeorology, 2012, 13, 1493–1506. http://dx.doi.org/10.1175/JHM-D-11-056.1. Armstrong, R. L., Chang, A., Rango, A., & Josberger, E. Snow depths and grain-size relationships with relevance for passive microwave studies. Annals of Glaciology, 1993, 17, 171–176. Boyarskii, D. A., & Etkin, V. S. Two flow model of wet snow microwave emissivity. IEEE International Geoscience and Remote Sensing Symposium 1994, Pasadena, CA, 1994, 4, 2068–2070. Brown, R., & D. Robinson. Northern Hemisphere spring snow cover variability and change over 1922–2010 including an assessment of uncertainty, The Cryosphere, 2011, 5(1), 219-229, http://dx.doi.org/10.5194/tc-5-219-2011. Butt, M. J., & Kelly, R. E. J. Estimation of snow depth in the UK using the HUT snow emission model. International Journal of Remote Sensing, 2008, 29(14), 4249–4267. http://dx.doi.org/10.1080/01431160801891754. Butt, M. J. A comparative study of Chang and HUT models for UK snow depth retrieval. International Journal of Remote Sensing, 2009, 30(24), 6361–6379. http://dx.doi.org/10.1080/01431160902852804. Chang, A. T. C., Foster, J. L., & Hall, D. K. Nimbus-7 SMMR derived global snow cover parameters. Annals of Glaciology, 1987, 9, 39–44. Chang, W. , Tan, S. , Lemmetyinen, J., Tsang, L., Xu, X., & Yueh, S. H. Dense media radiative transfer applied to SnowScat and SnowSAR. IEEE J. Sel. Top. Appl. Earth Observations Remote Sensing, 2014, 7(9), 3811–3825. Clifford, D. Global estimates of snow water equivalent from passive microwave instruments: history, challenges and future developments. International Journal of Remote Sensing, 2010, 31(14), 3707– 3726.http://dx.doi.org/10.1080/01431161.2010.48348 2. 133 Cline, D., R. Armstrong, R. Davis, K. Elder, & G. Liston. CLPX LSOS snow pit measurements. Edited by M. Parsons and M.J. Brodzik. In CLPX-Ground: Snow measurements at the Local Scale Observation Site (LSOS), J. Hardy, J. Pomeroy, T. Link, D. Marks, D. Cline, K. Elder, R. Davis. 2003. Boulder, CO: NASA DAAC at the National Snow and Ice Data Center, 2003, Updated July 2004. Chen, C. T., Tsang, L., Guo, J., Chang, A. T. C., & Ding, K. H. Frequency dependence of scattering and extinction of dense media based on three-dimensional simulations of Maxwell's equations with applications to snow. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(8), 1844–1852. Cui, Y., Xiong, C., Lemmetyinen, J., Shi, J., Jiang, L., Peng, B., Li, H., Zhao, T., Ji, D., & Hu, T. Estimating snow water equivalent with backscattering at X and Ku band based on absorption loss. Remote Sensing, 2016, 8(6), 505–18. http://doi.org/10.3390/rs8060505. Davis, D. T., Chen, Z., Tsang, L., Hwang, J. N., & Chang, A. T. C. Retrieval of snow parameters by iterative inversion of a neural network. IEEE Transactions on Geoscience and Remote Sensing, 1993, 31, 842–852. http://dx.doi.org/10.1109/36.239907. Derksen, C., Toose, P., Rees, A., Wang, L., English, M., Walker, A., & Sturm, M. Development of a tundra-specific snow water equivalent retrieval algorithm for satellite passive microwave data. Remote Sensing of Environment, 2010, 114(8), 1699–1709. http://doi.org/10.1016/j.rse.2010.02.019. Derksen, C., Toose, P., Lemmetyinen, J., Pulliainen, J. T., Langlois, A., Rutter, N., & Fuller, M. C. Evaluation of passive microwave brightness temperature simulations and snow water equivalent retrievals through a winter season. Remote Sensing of Environment, 2012, 117, 236–248. Ding, K. H., Xu, X., & Tsang, L. Electromagnetic scattering by bicontinuous random microstructures with discrete permittivities. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(8), 3139–3151. Dobson, M. C., Ulaby, F. T., Hallikainen, M. T., & El-Rayes, M. A. Microwave dielectric behavior of wet soil-Part II: Dielectric mixing models. IEEE Transactions on Geoscience and Remote Sensing, 1985, (1), 35–46. http://dx.doi.org/10.1109/TGRS. 1985.289498. Du, J., Shi, J., & Xiong, C. A method to estimate snow water equivalent using multiangle X-band radar observations (pp. 3774–3776). Presented at the IEEE International Geoscience and Remote Sensing Symposium (IGARSS) 2010, Honolulu, HI. 2010a. http://doi.org/10.1109/IGARSS.2010.5651642. 134 Du, J., Shi, J., & Rott, H. Comparison between a multi-scattering and multi-layer snow scattering model and its parameterized snow backscattering model. Remote Sensing of Environment, 2010b, 114(5), 1089–1098. http://doi.org/10.1016/j.rse.2009.12.020. Durand, M., Kim, E. J., & Margulis, S. A. Quantifying uncertainty in modeling snow microwave radiance for a mountain snowpack at the point-scale, including stratigraphic effects,” IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(6), 1753– 1767. Durand, M., Kim, E. J., Margulis, S. A., & Molotch, N. P. A first-order characterization of errors from neglecting stratigraphy in forward and inverse passive microwave modeling of snow. IEEE Geoscience and Remote Sensing Letters, 2011, 8(4), 730–734. http://dx.doi.org/10.1109/LGRS.2011.2105243. Durand, M., & Liu, D. The need for prior information in characterizing snow water equivalent from microwave brightness temperatures. Remote Sensing of Environment, 2012, 126(C), 248–257. http://doi.org/10.1016/j.rse.2011.10.015. ESA. Report for Mission Selection: CoReH2O, ESA SP-1324/2 (3 volume series), European Space Agency (ESA), Noordwijk, The Netherlands, 2012. Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P.K. and Sokratov, S.A. The international classificationfor seasonal snow on the ground. IHP-VII Technical Documents in Hydrology N° 83 | IACS Contribution N° 1, UNESCO-IHP, Paris, 2009, 1–80. Flanner, M., K. Shell, M. Barlage, D. Perovich, & Tschudi, M. Radiative forcing and albedo feedback from the Northern Hemisphere cryosphere between 1979 and 2008, Nature Geoscience, 2011, 4(3), 151-155, doi:10.1038/NGEO1062. Frei, A., Tedesco, M., Lee, S., Foster, J., Hall, D. K., Kelly, R., & Robinson, D. A. A review of global satellite-derived snow products. Advances in Space Research, 2012, 50(8), 1007–1029. http://dx.doi.org/10.1016/j.asr.2011.12.021. Gan, T. Y., Kalinga, O., & Singh, P. Comparison of snow water equivalent retrieved from SSM/I passive microwave data using artificial neural network, projection pursuit and nonlinear regressions. Remote Sensing of Environment, 2009, 113(5), 919–927. http://doi.org/10.1016/j.rse.2009.01.004. Graf, T., Koike, T., Fujii, H., Brodzik, M., & Armstrong, R. L. CLPX-ground: ground based passive microwave radiometer (GBMR-7) data. Boulder, Colorado USA: NASA DAAC at the National Snow and Ice Data Center, 2003. 135 Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. Metropolis and MetropolisHasting algorithms. In Bayesian Data Analysis (2nd Edition) (Chapter 11.2). Boca Raton, FL: Chapman & Hall/CRC. 1995, 320–334. Grubbs, F. Proceduresfor detecting outlying observationsin samples. Technometrics, 2007, 11(1), 1–21. http://dx.doi.org/10.1080/00401706.1969.10490657. Hallikainen, M. T., Ulaby, F. T., & Van Deventer, T. E. Extinction behavior of dry snow in the 18-to 90-GHz range. IEEE Transactions on Geoscience and Remote Sensing, 1987, 25(6), 737–745. Hardy, J., Pomeroy, J., Link, T., Marks, D., Cline, D., Elder, K., & Davis, R. CLPXground: snow measurements at the Local Scale Observation Site (LSOS). Boulder, Colorado USA: NASA DAAC at the National Snow and Ice Data Center, 2003. IGOS. Integrated Global Observing Strategy Cryosphere Theme Report - For the monitoring of our environment from space and from earth (No. WMO/TD-No. 1405). World Meteorological Organization, Geneva, Switzerland, 2007. Jiang, L., Shi, J., Tjuatja, S., Dozier, J., Chen, K., & Zhang, L. A parameterized multiplescattering model for microwave emission from dry snow. Remote Sensing of Environment, 2007, 111(2-3), 357–366. http://dx.doi.org/10.1016/j.rse.2007.02.034. Jiang, L., Wang, P., Zhang, L., Yang, H., & Yang, J. Improvement of snow depth retrieval for FY3B-MWRI in China. Science China: Earth Sciences, 2014, 57(6): 12781292. Jordan, R. A one-dimensional temperature model for a snow cover technical documentation for SNTHERM.89. Cold Regions Research & Engineering Laboratory. Special Report 91-16. 1991, 43–44. Joseph, J. H., Wiscombe, W. J., & Weinman, J. A. The delta-Eddington approximation for radiative flux transfer. Journal of the Atmospheric Sciences, 1976, 33, 2452–2459. Kelly, R. E., Chang, A. T., Tsang, L., & Foster, J. L. A prototype AMSR-E global snow area and snow depth algorithm. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(2), 230–242. http://dx.doi.org/10.1109/TGRS.2003.809118. Kontu, A., & Pulliainen, J. T. Simulation of spaceborne microwave radiometer measurements of snow cover using in situ data and brightness temperature modeling. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(3), 1031–1044. 136 Lemmetyinen, J., Pulliainen, J. T., Rees, A., Kontu, A., Qiu, Y., & Derksen, C. Multiplelayer adaptation of HUT snow emission model: comparison with experimental data. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(7), 2781–2794. Lemmetyinen, J., Kontu, A., Pulliainen, J., Vehviläinen, J., Rautiainen, K., Wiesmann, A., Matzler, C., Werner, C., Rott, H., Nagler, T., Schneebeli, M., Proksch, M., Schüttemeyer, D., Kern, M., Davidson, M. The Nordic Snow Radar Experiment. Geoscientific Instrumentation, Methods and Data Systems. 2016, 5, 403-415. Lettenmaier, D. P., Alsdorf, D., Dozier, J., Huffman, G. J., Pan, M., & Wood, E. F. Inroads of remote sensing into hydrologic science during the WRR era, Water Resources Research, 2015, 51(9), 7309-7342, doi:10.1002/2015WR017616. Liang, D., Xu, X., Tsang, L., Andreadis, K. M., & Josberger, E. G. The effects of layers in dry snow on its passive microwave emissions using Dense Media Radiative Transfer theory based on the Quasicrystalline Approximation (QCA/DMRT). IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(11), 3663–3671. http://dx.doi.org/10.1109/TGRS.2008.922143. Lowe, H., & Picard, G. Microwave scattering coefficient of snow in MEMLS and DMRT-ML revisited: the relevance of sticky hard spheres and tomography-based estimates of stickiness. The Cryosphere Discuss., 2015, 9(2), 2495–2542. Matzler, C., Schanda, E., & Good, W. Towards the definition of optimum sensor specifications for microwave remote sensing of snow. IEEE Transactions on Geoscience and Remote Sensing, 1982, GE-20, 57−66. http://dx.doi.org/10.1109/TGRS.1982.4307521. Matzler, C. Applications of the interaction of microwaves with the natural snow cover. Remote Sensing Reviews, 1987, 2(2), 259–387. Matzler, C. Improved Born approximation for scattering of radiation in a granular medium. Journal of Applied Physics, 1998, 83(11), 6111–6117. Matzler, C. Autocorrelation functions of granular media with free arrangement of spheres, spherical shells or ellipsoids. Journal of Applied Physics, 1997, 81(3), 1509– 1517. Matzler, C., & Wiesmann, A. Extension of the microwave emission model of layered snowpacks to coarse-grained snow. Remote Sensing of Environment, 1999, 70(3), 317– 325. 137 Matzler, C. Mie Scattering with and without diffraction. Institute of Applied Physics, University of Bern, Research Report 2004-02, 2004a. Matzler, C. Notes on microwave radiation from snow samples and emission of layered snowpacks,” Institute of Applied Physics, University of Bern, Bern, Switzerland, IAP Research Report 96-09, 2004b. Matzler, C. Section 4.2: Comparison of emission models for covered surfaces. In Thermal Microwave Radiation: Applications for Remote Sensing, IET Electromagnetic Waves Series 52, London, UK: The Institution of Engineering and Technology, 2006, 227–240. Mironov, V. L., De Roo, R. D., & Savin, I. V. Temperature-dependable microwave dielectric model for an Arctic Soil. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(6), 2544–2556. http://doi.org/10.1109/TGRS.2010.2040034. Mocko, D., NASA/GSFC/HSL. NLDAS VIC Land Surface Model L4 Monthly Climatology 0.125 x 0.125 degree, version 002. Greenbelt, Maryland, USA: Goddard Earth Sciences Data and Information Services Center (GES DISC), 2014, Accessed June 2015. http://dx.doi.org/10.5067/B0L38R609B0R. Montpetit, B., Royer, A., Roy, A., Langlois, A., & Derksen, C. Snow microwave emission modeling of ice lenses within a Snowpack using the Microwave Emission Model for Layered Snowpacks. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(9), 4705–4717. Montpetit, B., Royer, A., Wigneron, J. P., & Chanzy, A. Evaluation of multi-frequency bare soil microwave reflectivity models. Remote Sensing of Environment, 2015, 162, 186–195. http://doi.org/10.1016/j.rse.2015.02.015. Nijseen, B., Schnur, R., & Lettenmaier, D. P. Global retrospective estimation of soil moisture using the Variable Infiltration Capacity Land Surface Model, 1980–93. Journal of Climate, 2001, 14, 1790–1808. http://dx.doi.org/10.1175/15200442(2001)014<1790:GREOSM>2.0.CO;2. Pan, J., Jiang, L., & Zhang, L. Comparison of the multi-layer HUT snow emission model with observations of wet snowpacks. Hydrological Processes, 2012, 28(3), 1071–1083. http://dx.doi.org/10.1002/hyp.9617. Pan, J., Durand, M., Vander-Jagt, B., & Liu, D. Application of the Markov Chain Monte Carlo algorithm for snow water equivalent retrieval from passive microwave measurements. Remote Sensing of Environment, 2017, 192: 150-165. 138 Picard, G., Brucker, L., Roy, A., Dupont, F., Fily, M., Royer, A., & Harlow, C. Simulation of the microwave emission of multi-layered snowpacks using the Dense Media Radiative transfer theory: the DMRT-ML model. Geoscientific Model Development, 2013, 6(4), 1061–1078. http://dx.doi.org/10.5194/gmd-6-1061-2013. Proksch, M., Matzler, C., Wiesmann, A., Lemmetyinen, J., Schwank, M., Lowe, H., & Schneebeli, M. MEMLS3&a: Microwave Emission Model of Layered Snowpacks adapted to include backscattering. Geoscientific Model Development, 2015, 8(8), 2611– 2626. http://doi.org/10.5194/gmd-8-2611-2015. Pulliainen, J. T., Grandell, J., & Hallikainen, M. T. HUT snow emission model and its applicability to snow water equivalent retrieval. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(3), 1378–1390. Pulliainen, J. T., & Hallikainen, M. T. Retrieval of regional snow water equivalent from space-borne passive microwave observations. Remote Sensing of Environment, 2001, 75(1), 76–85. http://dx.doi.org/10.1016/S0034-4257(00)00157-7. Pulliainen, J. T. Mapping of snow water equivalent and snow depth in boreal and subarctic zones by assimilating space-borne microwave radiometer data and ground-based observations. Remote Sensing of Environment, 2006, 101(2), 257–269. http://dx.doi.org/10.1016/j.rse.2006.01.002. Rodell, M. The SnowEx Campaign. NASA Goddard Space Flight Center, http://neptune.gsfc.nasa.gov/hsb/index.php?section=322, 2016. Roy, V., Goita, K., Royer, A., Walker, A. E., & Goodison, B. E. Snow water equivalent retrieval in a Canadian boreal environment from microwave measurements using the HUT snow emission model. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(9), 1850–1859. http://dx.doi.org/10.1109/TGRS. 2004.832245. Rott, H. Prospects of microwave remote sensing for snow hydrology. Hydrologie Applications of Space Technology. IAHS Pul., no, 160, 1986. Rott, H., Nagler, T., Voglmeier, K., Kern, M., Macelloni, G., Gai, M., Cortesi, U., Scheiber, R., Hajnsek, I., Pulliainen, J., Flach, D. Algorithm for retrieval of snow mass from Ku- and X-band radar backscatter measurements. Presented at the IEEE International Geoscience and Remote Sensing Symposium (IGARSS) 2012, Munich. 2012, 135-188. http://doi.org/10.1109/IGARSS.2012.6350911. Roy, A., Picard, G., Royer, A., Montpetit, B., Dupont, F., Langlois, A., Derksen, C., & Champollion, N. Brightness temperature simulations of the Canadian seasonal snowpack driven by measurements of the snow Specific Surface Area. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51, 9, 4692–4704. 139 Royer A., Roy, A., Montpetit, B., Saint-Jean-Rondeau, O., Picard, G., Brucker, L., & Langlois, A. Comparison of commonly-used microwave radiative transfer models for snow remote sensing. Remote Sensing of Environment, 2017, 190, 247– 259, http://dx.doi.org/10.1016/j.rse.2016.12.020. Rutter, N., Sandells, M., Derksen, C., Toose, P., Royer, A., Montpetit, B., Langlois, A., Lemmetyinen, J., & Pulliainen, J. T. Snow stratigraphic heterogeneity within groundbased passive microwave radiometer footprints: Implications for emission modeling. Journal of Geophysical Research - Earth Surface, 2014, 119(3), 550–565. Schneebeli M. & Sokratov, S. A. Tomography of temperature gradient metamorphism of snow and associated changes in heat conductivity. Hydrol. Process., 2004, 18(18), 3655– 3665 Shi, J., Dozier, J. Estimation of snow water equivalence using SIR-C/X-SAR, Part II: Inferring Snow Depth and Particle Size, 2000, 38(6), 2475–2488. Sturm, M., Holmgren, J., & Liston, G. E. A seasonal snow cover classification system for local to global applications. Journal of Climate, 1995, 8(5), 1261–1283. http://dx.doi.org/10.1175/1520-0442(1995)008<1261:ASSCCS>2.0.CO;2. Sturm, M., Taras, B., Liston, G. E., Derksen, C., Jonas, T., & Lea, J. Estimating snow water equivalent using snow depth data and cimate classes. Journal of Hydrometeorology, 2010, 11(6), 1380–1394. http://dx.doi.org/10.1175/2010JHM1202.1. Takala, M., Luojus, K., Pulliainen, J. T., Derksen, C., Lemmetyinen, J., Kärnä, J. P., Koskinen, J., & Bojkov, B. Estimating northern hemisphere snow water equivalent for climate research through assimilation of space-borne radiometer data and ground-based measurements. Remote Sensing of Environment, 2011, 115(12), 3517–3529. http://dx.doi.org/10.1016/j.rse.2011.08.014. Tan, S., Chang, W., & Tsang, L. Modeling both active and passive microwave remote sensing of snow using dense media radiative transfer (DMRT) theory with multiple scattering and backscattering enhancement. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(9), 4418–4430. http://doi.org/10.1109/jstars.2015.2469290. Tedesco, M., Kelly, R., Foster, J. L., & Chang, A. T.C. AMSR-E/Aqua Daily L3 global snow water equivalent EASE-Grids. Version 2. Boulder, Colorado USA: NASA National Snow and Ice Data Center Distributed Active Archive Center. 2004a. doi: 10.5067/AMSR-E/AE_DYSNO.002. 140 Tedesco, M., Pulliainen, J. T., Takala, M., Hallikainen, M. T., & Pampaloni, P. Artificial neural network-based techniques for the retrieval of SWE and snow depth from SSM/I data. Remote Sensing of Environment, 2004b, 90(1), 76–85. http://dx.doi.org/ 10.1016/j.rse.2003.12.002. Tedesco M. & Kim, E. J. Intercomparison of electromagnetic models for passive microwave remote sensing of snow. IEEE Transactions on Geoscience and Remote Sensing, 2006a, 44(10), 2654–2666. Tedesco, M., & Kim, E. J. Retrieval of dry-snow parameters from microwave radiometric data using a dense-medium model and genetic algorithms. IEEE Transactions on Geoscience and Remote Sensing, 2006b, 44(8), 2143–2151. http://dx.doi.org/10.1109/TGRS.2006.872087. Tedesco, M. NRT AMSR2 Daily L3 global snow water equivalent EASE-grids. Dataset available online from NASA LANCE AMSR2 at the GHRC DAAC Huntsville, Alabama, U.S.A. 2015. doi: http://dx.doi.org/10.5067/AMSR2/A2_DySno_NRT. Tsang, L., Ding, K. H., & Wen, B. Dense media radiative transfer theory for dense discrete random media with particles of multiple sizes and permittivities. Progress in Electromagnetic Research, 1992, 6(5), 181–225. Tsang, L., Chen, C. T., Chang, A. T. C., Guo, J., & Ding, K. H. Dense media radiative transfer theory based on quasicrystalline approximation with applications to passive microwave remote sensing of snow,” Radio Science, 2000, 35(3), 731–749. Tsang, L. & Kong, J. A. Multiple scattering theory for discrete scatterers. In Scattering of electromagnetic waves - Advanced topics (Chapter 5). New York: John Wiley & Sons. Inc. 2001, 197-241. Tsang, L., Pan, J., Liang, D., Li, Z., Cline, D. W., & Tan, Y. Modeling active microwave remote sensing of snow using Dense Media Radiative Transfer (DMRT) theory with multiple-scattering effects. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(4), 990–1004. Ulaby, F. T., & Long, D. G. Emission by snow-covered terrain. In Emission models and land observations, Microwave Radar and Radiometric Remote Sensing (Chapter 12-11). Ann Arbor, MI, USA: The University of Michigan Press. 2014, 586-595. Wang, H. , Pulliainen, J. T., & Hallikainen, M. T. Application of strong fluctuation theory to microwave emission from dry snow. Progress in Electromagnetics Research, 2000, 29, 39–55. 141 Wegmuller, U., & Matzler, C. Rough bare soil reflectivity model. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(3), 1391–1395. http://dx.doi.org/10.1109/36.763303 Wiesmann, A., Matzler, C., & Weise, T. Radiometric and structural measurements of snow samples. Radio Science, 1998, 33(2), 273–289. Wiesmann, A., & Matzler, C. Microwave emission model of layered snowpacks. Remote Sensing of Environment, 1999, 70(3), 307–316. Wiesmann, A., & Matzler, C. Technical Documentation and Program Listings for MEMLS 99.1, Microwave emission model of layered snowpacks. Institute of Applied Physics, University of Bern, IAP Report Nr. 99-4, 15.09.1999, 1999b. Xiong, C., Shi, J., & Lemmetyinen, J. Refinement of the X and Ku band dual-polarization scatterometer snow water equivalent retrieval algorithm. Presented at the IEEE International Geoscience and Remote Sensing Symposium (IGARSS) 2014, Cape Town, South Africa. 2014, 2419–2422. http://doi.org/10.1109/igarss.2014.6946960. Yueh, S. H., Dinardo, S. J., Akgiray, A., West, R., Cline, D. W., & Elder, K. Airborne Ku-Band polarimetric radar remote sensing of terrestrial snow cover. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(10), 3347–3364. http://doi.org/10.1109/TGRS.2009.2022945. 142 Appendix A: Solving A0 and B0 For 1-flux RT Equation Given boundary conditions as Tup(-d’+) and Tdn(0-), the constants, A0 and B0, can be solved by substituting them into the general solution, as: Tup (−d ') = T + γ 0 A0 exp(−γ d ') + 1 B exp(γ d ') γ0 0 Tdn (0) = T + A0 + B0 (A1-1) (A1-2) Then, A0 and B0 can be solved as: A0 = "T (0− ) − T $ 1 exp(γ d ') − "T (−d '+ ) − T $ # dn %γ # up % 0 1 exp(γ d ') − γ 0 exp(−γ d ') γ0 (A1-3) "T (−d '+ ) − T $ − "T (0− ) − T $γ exp(−γ d ') # up % # dn % 0 B0 = 1 exp(γ d ') − γ 0 exp(−γ d ') γ0 143 (A1-4) Appendix B: Solving RT Equations For Multi-Layer Snowpack This appendix shows how the RT solutions for homogenous snow case, i.e., equation (2.13)-(2.14) for MEMLS and (2.21)-(2.22) for the HUT model, could be used to derive the published equations for multiple snow layers. As shown in Figure 2, define the four fluxes, Aj, Bj, Cj, Dj as the downward and upward TB just above the lower boundary (at z = zj-1), and the downward and upward TB just below the upper boundary (at z = zj) of the j-th snow layer, respectively. First, the relationships between A, B, C and D at interfaces between different layers are: B j = s j−1 Aj + (1− s j−1 )D j−1 (A2-1) C j = (1− s j )Aj+1 + s j D j (A2-2) where sj is the interface reflectivity of the upper boundary of the j-th snow layer; sj-1 is the interface reflectivity of the lower boundary. Applying equation (2.21)-(2.22) of the 1-flux model gives two additional equations as: D j = B j l j + K 0T j (1− l j ) (A2-3) Aj = C j l j + K 0T j (1− l j ) (A2-4) with, l j = exp "#−(κ e − qκ s )d ' j $% (A2-5) 144 K0 = κa κ e − qκ s (A2-6) where, Tj is the physical temperature of the j-th snow layer. d’j is the slant thickness of the j-th snow layer, defined as dj’= dj secθ0, where dj is the vertical thickness of this layer, θ0 is the observing angle. Instituting (A2-1) and (A2-2) into (A2-3) and (A2-4) gives, D j = "# s j−1 Aj + (1− s j−1 )D j−1 $% l + K 0T j (1− l j ) (A2-7) Aj = "#(1− s j )Aj+1 + s j D j $% l + K 0T j (1− l j ) (A2-8) Substituting equation (A2-8) into (A2-7) can eliminate Aj, and gives the following expression for Dj as: "(1− s )l D + s (1− s )l 2 A % j−1 j j−1 j−1 j j j+1 $ ' Dj = 2 1− s j−1s j l j $# + (1+ s j−1l j )(1− l j )K 0T j '& 1 (A2-9) Instead, substituting equation (A2-7) into (A2-8) derives the expression for Aj, as, "(1− s )l A + (1− s )s l 2 D % 1 j j j+1 j−1 j j j−1 $ ' Aj = 2 1− s j s j−1l j $# + (1+ s j l j )(1− l j )K 0T j '& (A2-10) Based on the definitions of A and D, It can be found that equation (A2-9) and (A2-10) are identical to equation (6) and (8) of multiple-layer HUT model in Lemmetyinen et al. (2010). Similarly, applying equation (2.13)-(2.14) of the 2-flux model to the j-th snow layer could give: 145 1 B γ0 0 (A2-11) Aj = T j + A0 exp(−γ d ' j ) + B0 exp(γ d ' j ) (A2-12) D j = T j + γ 0 A0 + where, (B j − T j ) − (C j − T j ) A0 = γ 0 exp(−γ d ' j ) − B0 = − 1 exp(γ d ' j ) γ0 1 exp(γ d ' j ) γ0 (A2-13) (B j − T j ) − (C j − T j )γ 0 exp(−γ d ' j ) γ 0 exp(−γ d ' j ) − 1 exp(γ d ' j ) γ0 (A2-14) Substituting equation (A2-13) and (A2-14) into equation (A2-11) and (A2-12) gives, D j = T j + γ 0 A0 + 1 B γ0 0 ! $ 1 # γ 0 − + exp(−γ d ' j ) − exp(γ d ' j ) & γ0 & = T j #1− # & 1 γ 0 exp(−γ d ' j ) − exp(γ d ' j ) & # γ0 " % 1 γ0 + Bj 1 γ 0 exp(−γ d ' j ) − exp(γ d ' j ) γ0 γ0 − +Cj exp(−γ d ' j ) − exp(γ d ' j ) γ 0 exp(−γ d ' j ) − 1 exp(γ d ' j ) γ0 (A2-15) 146 Aj = T j + A0 exp(−γ d ' j ) + B0 exp(γ d ' j ) ! $ 1 # γ 0 − + exp(−γ d ' j ) − exp(γ d ' j ) & γ0 & = T j #1− # & 1 γ 0 exp(−γ d ' j ) − exp(γ d ' j ) & # γ0 " % 1 γ0 +Cj 1 γ 0 exp(−γ d ' j ) − exp(γ d ' j ) γ0 γ0 − exp(−γ d ' j ) − exp(γ d ' j ) + Bj γ 0 exp(−γ d ' j ) − 1 exp(γ d ' j ) γ0 (A2-16) Since equation (A2-15) and (A2-16) stand in general and they have several expressions in common, it is convenient to simplify them as, D j = e jT j + t j B j + rjC j (A2-17) Aj = e jT j + t jC j + rj B j (A2-18) where, 1 γ0 1− γ 02 tj = = t0 1 1− γ 02t02 γ 0 exp(−γ d ' j ) − exp(γ d ' j ) γ0 γ0 − rj = exp(−γ d ' j ) − exp(γ d ' j ) γ 0 exp(−γ d ' j ) − 1 exp(γ d ' j ) γ0 = (A2-19) 1− t02 γ0 1− γ 02t02 (A2-20) and, t0 = exp(−γ d ' j ) (A2-21) 147 The formulation of tj, rj and ej are solved rigorously, however, according to the roles they take in RT equation, they are actually “transmissivity”, “reflectivity” and “emissivity” of the j-th snow layer; t0 is called “one-way transmissivity through the snow layer” in Wiesmann & Matzler (1999a). Equation (2-17) to (2-21) are identical to the ones in MEMLS (Wiesmann & Matzler, 1999a). In multiple-layer HUT model (Lemmetyinen et al., 2010) and MEMLS (Wiesmann & Matzler, 1999a), given the equations that link the upward and downward TB of each snow layer to its adjacent layers, all the relationships and boundary conditions can be included in a big linear matrix. Finally, TB of all snow layers can be solved. Taking the upward TB just below surface can then easily calculate the TB observed by ground-based radiometer. More details of the entire process refer to the Appendix A of (Lemmetyinen et al., 2010) and the Documentation of MEMLS (Wiesmann & Matzler, 1999b). 148

1/--страниц