close

Вход

Забыли?

вход по аккаунту

?

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
Документ
Категория
Без категории
Просмотров
0
Размер файла
8 948 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа