JFY2015
Japan Aerospace Exploration Agency
Commissioned Research Agreement

Algorithm Theoretical Basis Document

GCOM-C1/SGLI Snow/Ice Products
Upgrades, testing and validation (GCOM-C1 4th RA 401)

JX-PSPC-419712

March 2016

Dr. Knut Stamnes (PI)Stevens Institute of Technology

SGLI Algorithm Description

1

CONTENTS

Contents
1 Introduction
1.1 Scope of the project . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Major objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Main goals and our work in JPY2015 . . . . . . . . . . . . . . . . . . . . . .
1.4 Overview of our work from JPY2013 to JPY2015 . . . . . . . . . . . . . . .

3
3
3
4
5

2 Theoretical Description of the Algorithms
2.1 C1: Cloud mask and surface classification . . . . . . . . . . . . . . . . . . .
2.1.1 Identifying snow/ice fields: the NDSI test . . . . . . . . . . . . . . .
2.1.2 Cloud tests over snow-covered land . . . . . . . . . . . . . . . . . . .
2.1.3 Cloud test over sea ice and ocean . . . . . . . . . . . . . . . . . . . .
2.1.4 Surface classification approach . . . . . . . . . . . . . . . . . . . . . .
2.1.5 Adaptation of C1 algorithm to SGLI sensor . . . . . . . . . . . . . .
2.2 C2: snow, sea ice and atmospheric parameter retrievals . . . . . . . . . . . .
2.2.1 Forward radiative transfer model . . . . . . . . . . . . . . . . . . . .
2.2.2 Neural Network (NN) Training . . . . . . . . . . . . . . . . . . . . .
2.2.3 Nonlinear optimal estimation . . . . . . . . . . . . . . . . . . . . . .
2.2.4 Snow and sea ice surface albedo retrieval . . . . . . . . . . . . . . . .
2.3 C2: Surface temperature retrieval . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Algorithm description . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Aerosol over snow retrieval by use of polarization channels . . . . . . . . . .

6
6
6
7
10
11
13
19
20
25
26
27
28
29
30
42

3 Algorithm Application and Validation
3.1 Application and validation of the C1 algorithm . . . . . . . . . . . . . . . . .
3.1.1 Image-based tests over Greenland . . . . . . . . . . . . . . . . . . . .
3.1.2 Direct comparison with CALIOP using MODIS Aqua data . . . . . .
3.1.3 Cross comparison of C1 using MODIS Terra and Aqua data . . . . .
3.1.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Application and validation of the C2 algorithm . . . . . . . . . . . . . . . . .
3.2.1 Comparison with GLI algorithms and field measurements . . . . . . .
3.2.2 Application to MODIS data . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 Snow albedo and temperature – comparisons with field measurement
data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44
44
46
47
49
50
52
52
55

References

63

SGLI Algorithm Description

2

60
61

1

Introduction

1.1

Scope of the project

To create GCOM-C1/SGLI snow/ice products we have developed algorithms for cloud mask
and surface classification, snow grain size retrieval, and surface temperature retrievals. Thesealgorithms are based on the work we did for ADEOS-II/GLI (Stamnes et al., 2007; Aoki etal., 20017; Hori et al., 2007), but we used the new SGLI channels, upgraded the surface classification algorithm, and improved the accuracy of the retrieved snow grain size and aerosol
properties as well as the surface temperature. These improvements have been accomplishedby using a linearized radiative transfer model (RTM) for the coupled atmosphere-snow system in conjunction with a 2-layer snow model to allow for simultaneous retrieval of snow andaerosol properties, and by using the response functions of the SGLI IR channels and updated
sounding data for the temperature retrieval. To test the performance of these algorithms
and thereby obtain further improvement, we proposed to (i) upgrade the quality and accuracy of our cloud mask algorithm over snow/ice cover area, (ii) upgrade and test our surface
classification algorithm over mixed snow/forest areas, (iii) develop and test algorithms for
direct retrieval of snow/sea ice albedo, (iv) explore the potential for improved snow retrieval
by assuming non-spherical snow grains, (v) use the direct albedo retrieval method to assess
the quality of albedo estimates derived from inferred snow parameters, and (vi) explore themerits of using polarization channels to improve the aerosol retrieval to the extent available
time and resources permit.
Figure 1 and Table 1 provide a summary of our products and the research schedule foreach of our product, respectively.

1.2

Major objectives

The major objectives of our work are to:
1. Develop an improved cloud mask algorithm over snow/ice covered areas that is optimized for high latitude and high surface elevation areas.
2. Develop an improved surface classification algorithm with improved snow/sea ice,
snow/forest classification.
3. Develop and test algorithms for direct retrieval of snow/sea ice albedo.
4. Explore the potential for improved snow retrieval by assuming non-spherical snow
grains.
5. Use the direct albedo retrieval method to assess the quality of albedo estimates derived
from infrared snow parameters.
6. Explore the merits of using polarization channels to improve the aerosol retrieval tothe extent available time and resources permit.

SGLI Algorithm Description3

1.3 Main goals and our work in JPY2015

L1Bdata!

1km/500m

1km/250m

L1B1 VN

L1B1 SWI

L1B1 TI

L1B1 P!"#$%!&'!"#$%&'%('")(*'+,$-"($)%!&'!"./-/%(,0")(*'+,$-!

1km

Daytime data!
Snow, Sea ice & Cloud area classification
1. Snow and Sea Ice Distribution!"#$%&'($"&)$*+,$+-"..*/+"0'&!
Snow & Sea ice classification
3. Snow and sea ice classification

2. Sea Ice Distribution in the Sea of Okhotsk!1#$234.*+"-$56'5,674$6,76*,8"-.!
Snow physical param. 1

RTM + Neural Network + Non-linear optimal estimation !

4. Snow surface temp.

Snow phys. param. 2

Snow phys. param. 3

5. Shallow Grainsize (0.86um)

6. Sub-sfc. Grain size
7. Surface grain size
8. Impurities

Atm. phys. param.
10. Aerosol optical
thickness

Forward RTM!
a. Snow metamorphism
potential map

b. Map of melt
onset date

Snow phys. param. 4
9. Spectral albedo

c. Initial parameter forland-surface model

Figure 1: Schematic illustration of our products and the flow of the retrieval process. Blueand red boxes are the standard and research products, respectively.

1.3

Main goals and our work in JPY2015

The main goals and our work in JPY2015 can be summarized as follows:• To test and validate our new C1 cloud mask algorithm for snow covered land areas;(see Section 2.1.2);• To test and improve the C1 surface classification algorithm over Arctic ocean for sea-ice
area by using MODIS images (see Section 2.1.3);• To test and validate the implementation of C1 cloud mask algorithm to SGLI sensor(see Section 2.1.5);• To implement Voronoi particle snow model into the SGLI snow retrieval code (seeSection 2.2.1 about snow model and its phase function);• To compare the retrieved snow grain size from both spherical and non-spherical particle
snow models using MODIS data (see Section 3.2.2);• To update snow/ice surface temperature retrieval based on the new SGLI channels
response function (see Section 2.3).

SGLI Algorithm Description4

1.4 Overview of our work from JPY2013 to JPY2015• To explore the potential of using SGLI polarization channels to improve aerosol retrieval over snow (See Section 2.4).
The research schedule for these tasks are listed in Table 1.
Table 1: Our products for snow and ice classification shown in Fig. 1.
No.
1
2
3
4
5
6

1.4

item
Upgrade cloud mask over snow/ice covered areas
Upgrade surface classification
Develop sea ice surface direct albedo
Use direct albedo to assess the quality of snow parameters retrieval
Explore non-spherical snow grains retrieval
Explore SGLI polarization channels to improve the aerosol retrieval

Research Schedule
Aug. 2013 - Mar. 2014
Aug. 2014 - Mar. 2015
Apr. 2014 - Mar. 2015
Aug. 2013 - Mar. 2016
Apr. 2014 - Mar. 2016
Apr. 2014 - Mar. 2016

Overview of our work from JPY2013 to JPY2015

In this 3-year project, we upgraded and validated our cloud mark algorithm and surface
classification algorithm; developed a sea-ice surface direct albedo algorithm; improved thesnow grain size retrieval algorithm by using a non-spherical Voronoi particle model insteadof the spherical particle model; and explored the potential for improving the aerosol retrieval
by using SGLI polarization channels. An overview of the completed tasks are as follows:• Upgrade of the cloud mask over snow/ice area. A new set of tables with dynamic
thresholds was established to improve the quality of the cloud detection. It includes
the consideration of solar/viewing geometry, atmospheric vertical structure, and surface elevation. This new algorithm has be applied to MODIS data, and validated bycomparison with CALIPSO overpasses. Although SGLI has fewer thermal channels
than MODIS, the new algorithm performs better than the MODIS MYD35 product[Chen et al., 2014].
• Upgrade of the cloud mask over land, desert, and snow mixed with vegetation areas.
Based on the CALIPSO and MODIS products, we generate a neural network table
to improve the quality of the cloud mask over such areas. This unique method is
expected to imrpove the cloud mask, specially in the winter season and over desert
areas. Comparisons with MODIS and CALIPSO show a significant improvement over
our previous algorithm. The performance in the winter season is much better thanthat of the MODIS products.
• The sea-ice coverage detection has been updated too, but further validation is needed.
We are searching for measurements that can be used for this propose.
• A sea-ice model has been implemented in our radiative transfer model. Based on it, wedeveloped a new algorithm to retrieve the direct sea-ice surface albedo. This productis similar to our direct snow albedo algorithm. The albedo will be directly retrieved
SGLI Algorithm Description5

from the satellite radiances. Applications to real data retrieval and validation will beneeded in the future.
• Comparisons between snow direct albedo and indirect albedo, can be used to assess
the quality of retrieved snow parameters. Since the quality of the cloud mask and thesurface classification has been greatly improved, we now find the direct and indirect
albedo values to be very close. This method can remove pixels with large bias in thesnow retrieval.
• In order to verify the snow model used in snow grain size retrieval, we compared themodeled snow BRDF with measurements obtained by the NASA Cloud Absorption
Radiometer (CAR) deployed by aircraft over Elson Lagoon, Barrow, Alaska, on April
7, 2008. These comparisons show that non-spherical Voronoi particles provides better
agreement between modeled and measured snow BRD than a spheric particle model.
The current Voronoi particle model still leaves room for improvement, which will be
pursued in our work over the next three years. The new non-spherical Voronoi particle
snow model has been implemented into the snow grain size retrieval code. Preliminary
result are presented in this report.
• We did a first sensitivity study using the polarization channels. The results showthat the polarized reflectance provides an opportunity to distinguish aerosol optical
thickness from snow impurity. The design of the SGLI polarized sensor helps makethis distinction possible, because the scattering angles are in the forward direction.
We have completed the 3-year RA4 project. The cloud mask algorithm C1 and the snow
retrieval algorithm C2 have been greatly improved and validated. For the next 3 years of
the RA6 project, we will focus on further validation of these two algorithms using MODIS
data and the new SGLI data. Also, we will improve the non-spherical particle snow modeland continue exploring the merit of using SGLI polarization channels and multiple-angle
measurements to improve aerosol retrievals over snow.

2

Theoretical Description of the Algorithms

2.1

C1: Cloud mask and surface classification

2.1.1

Identifying snow/ice fields: the NDSI test

The processing flow of the C1 algorithm starts with a Normalized Difference Snow Index(NDSI) test. The NDSI method has a long history and has been extensively utilized forcloud screening and detection of snow/ice covered areas [Hall et al., 1995, 1996, 2001]. The
NDSI is defined as a normalized difference between the reflectance in two bands:
NDSI =RVIS − RSWIRRVIS + RSWIR(1)where RVIS is the observed reflectance of a channel in the visible spectral range and RSWIR
is the observed reflectance of a channel in the shortwave infrared spectral range. For the
SGLI Algorithm Description

62.1 C1: Cloud mask and surface classification

MODIS instrument RVIS is selected to be band 4 (0.55 µm) and RSWIR to be band 6 (1.64µm, Terra) or band 7 (2.13 µm, Aqua). By choosing an appropriate threshold (NDSI ≥ 0.4
as proposed by Hall et al. [1995]), the NDSI method can effectively identify snow/ice covered
areas from other surface types as shown by Ackerman et al. [1998] and Salomonson & Appel[2004, 2006]. It can also be used to separate clouds (especially water clouds) from snow/ice
surfaces. In the C1 algorithm we apply the the NDSI method of Hall et al. [2001] because
we have found it to be very effective in identifying possible snow/ice covered areas including
partially snow-covered areas (snow mixed with vegetation). However, we have also found,
as discussed by Hutchison et al. [2013], that the NDSI method alone cannot completely
separate clouds from snow, because some clouds, especially thin ice clouds and overlapping
ice and water clouds, have NDSI signatures very similar to those of snow and ice. Therefore,in our C1 algorithm the NDSI method is used as a first step to identify possible snow/ice
covered areas for which more detailed cloud mask procedures must be applied for correct
cloud identification. A flowchart of the algorithm is provided in Figure 2.

Satellite	
  Radiance	
  Data	
  

NDSI	
  
Test	
  

not snow/ice

snow/ice possible•••

Solar/viewing	
  geometry	
  Atmospheric	
  proﬁle	  Surface	  eleva@on	  SWIR	
  Dynamic	  Threshold	
  Tests	
  Cloud	
  Tests	
  	
  
Over	
  
Land/Ocean	  CLOUD	  Snow/ice	  Classiﬁca@on	
  

SNOW	
  

SNOW	
  

SEA	
  ICE	  Over	
  vegeta/on	
  

Over	
  land	  With	
  snow	
  cover	
  

LAND	
  

OCEAN	  SEA	  ICE	
  Figure 2: Flow chart of the C1 cloud mask algorithm over snow/ice surfaces.

2.1.2

Cloud tests over snow-covered land

The use of SWIR tests to separate clouds from snow has been discussed in some detail byMiller & Lee [2005], who used a combination of VIS, SWIR and 1.38 µm thresholds with asimple solar zenith angle correction. Trepte et al. [2002] used the MODIS 1.64 µm reflectanceas well as the reflectance ratio between the 1.64 µm and the 0.65 µm channels to distinguishclouds from snow in polar regions. We tried to use a single SWIR threshold (2.13 µm) firstto distinguish clouds from snow surfaces and found that most liquid water clouds and thickice clouds have high 2.13 µm reflectance (R2.13 ) because of their small particle size or high
SGLI Algorithm Description

72.1 C1: Cloud mask and surface classification

location in the atmosphere, which makes thresholding possible. However, we also found thata fixed threshold of R2.13 is not appropriate because of the anisotropy of snow reflectance inthe 2.13 µm channel. At low solar illumination angles (i.e. solar zenith angles ≥ 70◦ ) the
threshold value should be higher than at nadir to avoid mis-classification of snow as clouds.
Furthermore, we found seasonal and surface elevation dependence of the threshold value
mainly due to different atmospheric absorption (primarily water vapor absorption) of thereflected 2.13 µm signal under different surface elevation and atmospheric conditions. For
example in Greenland, snow surfaces generally look “brighter” in Spring than in Summer
and Autumn, implying that higher thresholds must be applied. Also, we should use ahigher threshold value at Automated Weather Station (AWS) Summit (elevation: 3,254 m)
than at Swiss camp (elevation: 1,149 m) since the snow surface at Summit would generally
have higher 2.13 µm reflectance than the same snow surface at Swiss camp under the same
solar/viewing geometry. In order to accommodate these many factors, we use the following
strategy: The threshold is set to be the ‘maximum possible snow reflectance’ defined as thesimulated R2.13 reflectance of ‘very fine grained’ snow for a given elevation and atmosphericcondition. We chose to model reflectance of snow instead of clouds for the cloud mask
development because there could be many types of clouds over snow making it difficult todecide what kind of R2.13 reflectance a cloud should have. However, if we can establish “themaximum R2.13 value that a natural snow surface can have” then the situation becomes moreclear. We finally chose R2.13 reflectance of snow with effective grain size reff ≈ 25 µm as thethreshold (ThR2.13 ) since this grain size is smaller than most of the snow grains that exist
in nature. By applying this threshold, water clouds and thick ice clouds with high 2.13 µm
reflectance can be distinguished without mis-classifying most snow pixels as clouds. Lookup
tables that cover all possible solar/viewing geometries applicable for daytime operation of
the SGLI were prepared so that for each satellite pixel an appropriate threshold value canbe dynamically determined consistent with the actual solar/viewing geometry.
We further explored the dynamic threshold method by using two SWIR channels (i.e.
MODIS channels at 1.64 µm and 2.13 µm for Terra, and at 1.24 µm and 2.13 µm for
Aqua). Figure 3 is a radiative transfer simulation of MODIS channel 6 (1.64 µm) andchannel 7 (2.13 µm) reflectances of water/ice and mixed phase clouds over snow surfaces
with different snow grain sizes under a common solar/viewing geometry (solar zenith angle(SZA) = 60◦ , viewing zenith angle (VZA) = 20◦ , relative azimuth angle (RAZ) = 112◦ ) in
the polar regions. In the simulations we used the Sub-Arctic Summer atmospheric profile;the snow surfaces are assumed to consist of snow with effective grain sizes of 50, 100 and
400 µm; a liquid water cloud, located at 2 km above the snow surface, assumed to have aneffective particle size of 15 µm; an ice cloud, located at 10 km altitude, assumed to havean effective particle size of 42 µm. The inherent optical properties (IOPs) for water clouds
were obtained from the parameterization of Hu & Stamnes [1993] and ice cloud IOPs weretaken to be that of “rough aggregate” particles in the parameterization by Key et al. [2002].
The optical depth for water/ice clouds at 0.645 µm ranged between 0.01 and 5. Additional
calculations of the snow cases (grain radius range: 15 - 2000 µm) in the absence of clouds
were also performed. All the calculations were performed using an extended version of the
DISORT [Stamnes et al., 1988a] radiative transfer model with a Earth curvature correction[Dahlback & Stamnes, 1991]. In Fig. 3 reflectances of clouds over snow cases are represented
SGLI Algorithm Description

82.1 C1: Cloud mask and surface classificationby colored symbols: green for ice clouds, cyan for water clouds and blue for mixed ice andwater clouds. Different symbols indicate different cloud optical depths. For example, thethree green ‘crosses’ represent the reflectance of ice clouds with optical depth τi = 0.1 oversnow with reff = 400, 100, and 50 µm, while the green ‘circles’ are for τi = 0.5 cases. Greendashed lines indicate ice cloud cases with different optical depths over an underlying snow
surface with the same grain size. Reflectances of snow in the absence of cloud cover are
represented by the black dots.
0.25

snow surface with different grain sizes (τw = 0, τi = 0)
model suggested thresholdwater clouds (τw = 0.01)
water clouds (τ = 0.1)

ice clouds (τi = 0.01) over water clouds (τw = 0.01)

0.2

ice clouds (τi = 0.1) over water clouds (τw = 0.01)
ice clouds (τi = 0.5) over water clouds (τw = 0.01)

r

snow

0.15

=

ice clouds (τi = 1.0) over water clouds (τw = 0.01)= 0 µm

ice clouds (τL = 5.0) over water clouds (τw = 0.01)

rsnow = 15 µm

m
0µ


ice clouds (τi = 0.1)
ice clouds (τi = 0.5)w

r snoice clouds (τi = 1.0)
ice clouds (τi = 5.0)µm

rsnow = 20 µm

ow =

40

0

0.1

rsnow = 25 µmr

sn

MODIS channel 7 (2.13 µm) reflectancew

water clouds (τw = 0.5)

rsnow = µm

0.05
rsnow = 0 µm
rsnow = µm

00200 µm < rsnow < 2000 µm

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

MODIS channel 6 (1.64 µm) reflectance

Figure 3: Simulated R2.13 and R1.64 of clouds over snow. rsnow is the effective grain radius of
the snow; three ‘background’ snow grain radii (50, 100 and 400 µm) are used when simulating
the reflectances of cloud over snow, τw and τi are water and ice cloud optical depth at 0.645µm respectively. Sub-Arctic Summer atmospheric profile is used in the simulation and thesolar/viewing geometry adopted in this figure is: solar zenith angle (SZA) = 60◦ , viewing
zenith angle (VZA) = 20◦ , relative azimuth angle (RAZ) = 112◦ .
From the simulated results we can infer:• The reflectance from water clouds in the SWIR channels increases drastically as their
optical depth increases and it saturates at a higher level than that of snow. More
importantly, the reflectance from water clouds depends little on surface properties(except for the extremely thin (τi = 0.01) cases), which might indicate that a constantR2.13 or R1.64 threshold would be good enough to separate water clouds from underlying
snow surfaces.
SGLI Algorithm Description

92.1 C1: Cloud mask and surface classification• The reflectances of mixed ice and liquid water clouds are largely dominated by thecontribution from water droplets (even with optical depth 0.01) implying that theymight also be distinguished from snow in a way similar to that for water clouds.
• For ice clouds the situation is more complicated. We can see that very thin ice clouds(like the green ‘crosses’ with τi = 0.1) are very difficult to distinguish from snow since
their reflectances closely resemble those of snow with smaller grain sizes. We also findthat the reflectance for thin to moderately thick ice clouds (τi = 0.1 → 1) depends
strongly on the reflectance contribution of the snow surfaces underneath. The green‘circles’ in Figure 3 show that ice clouds with τ ≈ 0.5 or thinner might be separated
from the cloud free cases by a proper thresholding method.
An illustration of the ‘Model Suggested Threshold’ (MST) for the particular solar/viewing
geometry adopted in Fig. 3 is shown as the red line. The threshold consists of two parts:
the horizontal part is the aforementioned R2.13 threshold (ThR2.13 ), which is calculated fromthe reflectance of ‘very fine grained’ snow (reff ≈ 25 µm) and is designed to pick up thickice clouds and liquid water clouds. The curved part (blue line) is a quadratic fit to the 1.64
µm and 2.13 µm reflectances of snow, designed to pick up thinner ice clouds. For a possible
2
snow pixel it is calculated as ThR2.13 1.64 = c0 + c1 R1.64 + c2 R1.64, where R1.64 is the 1.64 µm
reflectance of that pixel, the symbol  indicates the dependence of this threshold on theR1.64 reflectance. The fitting coefficients c0 , c1 and c2 depend on the solar/viewing geometry.
Figure 3 shows that cloudy cases are ‘above’ the MST line in the R1.64 - R2.13 space which
implies that a snow pixel with 2.13 µm reflectance greater than the corresponding MST
value is suspected to be cloud contaminated. The ‘clear confidence’ indices are estimatedlinearly using the two parts of the MST (ThR2.13 and ThR2.13 1.64 respectively), ranging
between 0% clear (‘confident cloudy’) and 100% clear (‘confident clear’) and the final clear
confidence level is taken to be the minimum value of the two calculated clear confidence
indices. The thresholds and corresponding confidence levels are summarized in Table 2. An
artificial offset f is added to the curved part of the threshold to accommodate errors fromthe model. Currently f is set to be 0.0016 for the test over the Greenland region. The
simulation using the 1.24 µm and 2.13 µm channels of MODIS Aqua yielded similar results,
except that a four degree polynomial was needed to fit the curved part of the suggested
threshold. So the corresponding formula for the curved part of the threshold would become:2
3
4
ThR2.13 1.24 = c0 + c1 R1.24 + c2 R1.24+ c3 R1.24+ c4 R1.24. A multi-linear interpolation scheme isused to dynamically interpolate the suggested threshold line according to the solar/viewing
geometry of each satellite pixel. The threshold depends on season, location and surface
elevation. Currently, three model atmospheric conditions: Sub-Arctic Summer, Sub-Arctic
Winter and Mid-latitude Winter are tested and a model-based scaling method (please see
Chen et al. [2014] for details) is applied to deal with potential issues associated with thesurface elevation.
2.1.3

Cloud test over sea ice and ocean

The above SWIR dynamic threshold tests are all designed for snow covered land areas. Overocean areas we use the ratio of R2.13 /R0.65 instead of R2.13 since for some clouds over ocean
SGLI Algorithm Description10

2.1 C1: Cloud mask and surface classification

Table 2: SWIR test thresholds for clear confidence level determination
Threshold part
Horizontal Part
Curved Part

Threshold value and clear confidence
Confident clear
Confident cloudy
R2.13 ≤ ThR2.13
R2.13 ≥ ThR2.13 + 0.05
R2.13 ≤ ThR2.13 1.64 + f R2.13 ≥ ThR2.13 1.64 + f + 0.05

the 2.13 µm signal would be too low to pass the R2.13 threshold, and hence the algorithm will
mis-identify those cloudy pixels as sea ice pixels. However, the ratio R2.13 /R0.65 would still
be high enough to allow a proper thresholding method to be applied. Again, this thresholdis dynamically calculated like in the R2.13 method.
2.1.4

Surface classification approach

Another feature of the C1 algorithm is that it provides a detailed classification of the ice
type and snow coverage over sea ice. As discussed by Stamnes et al. [2011], different ice
types have very different ice albedo due to different amounts of scattering inclusions and ice
impurity concentrations, and small amounts of snow on sea ice will significantly change thealbedo of the ice. Thus, a detailed surface classification of sea ice is necessary to obtain arobust surface retrieval algorithm. We applied the Normalized Difference Ice Index (NDII)
method to achieve this goal in our GLI products [Stamnes et al., 2007a]. Based on our
experience with the GLI algorithm, we are developing a new surface classification algorithm
in C1 using our comprehensive radiative transfer model for the coupled atmosphere-snowice-ocean system [Stamnes et al., 2011].
Top of atmosphere sea ice and open water albedo and relative position of SGLI channels
1

open water
new ice (1cm)nilas (3cm)
nilas (5cm)grey ice (15cm)
first year sea ice (30cm)
SGLI channel VL05 (0.53 µ m)
SGLI channel VL07 (0.67 µ m)
SGLI channel VL10 (0.87 µ m)
SGLI channel SW01 (1.05 µ m)
SGLI channel SW03 (1.63 µ m)

0.9

0.8

Spectral albedo

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

400

600

800

1000

1200
1400
Wavelength [nm]

1600

1800

2000

2200

Figure 4: Simulated ice and open water spectral albedo and the relative position of SGLI
channels.
Figure 4 shows the spectral albedo of different types of sea ice as well as open water. It
can be seen that SGLI channel VL05 (0.53 µm) and VL07 (0.67 µm) are very sensitive tothe thickness of sea ice and thus would be the optimum channels for sea ice classification and
SGLI Algorithm Description

11

2.1 C1: Cloud mask and surface classification

ice/open water separation. Another important task of the C1 algorithm is to separate snowcovered sea ice pixels from bare sea ice pixels since it is very important to provide correct
pixel information to the C2 algorithm for sea ice/snow property retrieval. We simulated
the spectral albedo of snow-covered and bare sea ice as shown in Fig. 5. It can be seenthat snow and ice have very different albedo values in the NIR and SWIR regions, so SGLI
channel VL10 (0.87 µm) and SW01 (1.05 µm) can provide valuable information to separatesnow covered sea ice from bare sea ice. Due to the strong anisotropy of sea ice and snow
reflectance in the NIR and SWIR region, dynamic thresholds that depend on solar/viewing
geometry will be applied for this purpose.
Top of atmosphere sea ice albedo and relative position of SGLI channels
1

first year sea ice (30 cm thickness)
first year sea ice + snow (100 µ m grain size)
first year sea ice + snow (500 µ m grain size)
first year sea ice + snow (2000 µ m grain size)
SGLI channel VL05 (0.53 µ m)
SGLI channel VL07 (0.67 µ m)
SGLI channel VL10 (0.87 µ m)
SGLI channel SW01 (1.05 µ m)
SGLI channel SW03 (1.63 µ m)

0.9

0.8

Spectral albedo

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0

400

600

800

1000

1200
1400
Wavelength [nm]

1600

1800

2000

2200

Figure 5: Simulated spectral albedo of snow covered/bare first year sea ice.

SGLI Algorithm Description

122.1 C1: Cloud mask and surface classification

2.1.5

Adaptation of C1 algorithm to SGLI sensor

After our successful tests on MODIS data as described in Chen et al. [2014], we adapted ouralgorithm to use the channel specifications of GCOM-C1/SGLI and tested its performanceusing the SGLI test data. To distinguish clouds over snow covered area, we proposed to use
SGLI channel SW03 (1.63 µm) and SW04 (2.20 µm) since they are the closest equivalents to
MODIS channel 6 (1.64 µm) and channel 7 (2.13 µm). Following the same methodology as
in our previous work, we simulated various cloud/cloud-free cases with different underlying
snow conditions. The result, however, was very different from what we got for the MODIS
channels. Figure 6 shows the simulation of pure ice clouds with different optical depths (greendots) over different snow cases (black dots) using MODIS (left) and SGLI (right) channel
specifications. One can observe that the previously distinguishable cloudy/cloud-free cases
under MODIS channel specifications are no longer distinguishable by SGLI as the cloudy andcloud-free cases are closely packed together. This result was unexpected since we thought,
as described in Chen et al. [2014], that the combination of SGLI channel SW03 and SW04
should be an optimum choice for cloud screening over snow covered areas. We seriously
investigated this issue in order to find the reason behind it since it can, if unresolved, render
our SWIR cloud masking scheme inapplicable to the SGLI sensor.
We first noticed different profiles of gaseous absorption in MODIS and SGLI. Figure 7
0.35

snow surface with different grain sizes
pure ice clouds (τi = 0.1,0.5,1.0,2.0,5.0)

SGLI channel SW4 (2.20 µm) reflectance

MODIS channel 07 (2.13 µm) reflectance

0.25

0.2

0.15

0.1

0.05

0

snow surface with different grain sizes
pure ice clouds (τi = 0.1,0.5,1.0,2.0,5.0)

0.3

0.3

0.25

0.2

0.15

0.1

0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

MODIS channel 06 (1.64 µm) reflectance

0

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

SGLI channel SW3 (1.63 µm) reflectance

Figure 6: Simulated SWIR reflectance of snow (black dots) and ice clouds over snow (greendots). The ice cloud is assumed to be located at 10 km altitude with an effective particle sizeof 42 µm and a “background” snow effective grain sizes of 50, 100 and 400 µm. Cloud-freesnow cases are assumed to have effective grain size between 15 to 2000 µm. Other settings ofthe simulation are identical to the simulation in our previous work [Chen et al., 2014]. Leftpanel: For MODIS sensor; Right panel: For SGLI sensor, note the cloudy cases (greendots) are no longer separable from snow cases (black dots) in the simulation for SGLI.

SGLI Algorithm Description

13

2.1 C1: Cloud mask and surface classificationshows significant change in layer absorption optical depth between the two SWIR channels
from MODIS to SGLI. However, as calculated that this change in gas absorption profile
will only lead to a 4 to 7% difference in the simulation results, which will not be enough to
explain the difference between the results for MODIS and SGLI as shown in Figure 6. So
there must be some other reason for this difference. Next, we focused on the IOPs of clouds
and snow that we used in the simulations. Figure 8 shows the single scattering albedo andasymmetry factor of snow particles in the SWIR spectral range. We can notice that:• Snow particles have higher average single-scattering albedo in MODIS channel 6(1.64µm) compared to that of MODIS channel 7 (2.13µm), which means that snow ismore absorptive in MODIS channel 7 compared to MODIS channel 6. As the grain
size becomes smaller (from 400µm to 100µm) the difference in single-scattering albedo
between MODIS channel 6 and channel 7 becomes larger. However, for the SGLI channels the differences in single-scattering albedo between two SWIR channels are muchsmaller compared to their MODIS counterparts.
• For the MODIS channels, the average asymmetry factor of snow particles is significantly
higher in channel 7 than in channel 6. This difference implies that the phase functionof snow will be more asymmetric in channel 7 (scatters more in the forward direction)compared to channel 6. For SGLI, however, the difference in asymmetry factor is also
small between two SWIR channels.
The difference in snow IOPs suggests that we should look into the refractive index of ice
in SWIR region. Figure 9 shows the refractive index of ice in the 0.9 to 2.4 µm wavelength
range. One can see that while the real part of refractive index decreases smoothly from about
1.3 to 1.25, the imaginary part, which represents the absorption property of ice, has a peakat around 2µm. The close position of MODIS channel 7 to this peak results in significantly
more absorption by ice particles in this channel compared to all the other SWIR channels
in MODIS and SGLI (MODIS channels 5, 6, SGLI channels SW01, SW03 and SW04),
which explains the difference of snow IOPs we found earlier. Moreover, since cloud particles(only ice cloud is discussed here since they are hard to detect compared to water clouds)
are generally smaller than snow particles, the difference in IOPs between an “absorptive”
channel (like MODIS channel 7) and a “non-absorptive” channel (like MODIS channel 6 or
SGLI channel SW03) will be more significant. Figure 10 shows plots of the single-scattering
albedo in the two sensor’s SWIR channels with green colored symbols indicating those for
MODIS (channel 6 and 7) and blue colored symbols those for SGLI (channel SW03 andSW04). It can be seen that for the MODIS channels the symbols indicating cloud particles
are distinct from the line for snow particles with different grain sizes whereas this behavior is
not true for SGLI. Hence the single-scattering albedo of ice cloud particles in the two SWIR
channels resemble those of snow particles with small grain size for SGLI, which explains thesimulation results for the SGLI sensor (Fig. 6 right panel) where the reflected signals fromice clouds are converging to the line indicating cloud-free snow cases as the cloud opticaldepth increases from 0.1 to 5.
Knowing that the difference in ice absorption is the main reason that we see different resultsin our simulations, we now try to follow the strategy of selecting one “absorptive” channel
SGLI Algorithm Description

142.1 C1: Cloud mask and surface classification

Total absorption optical depth by atmospheric gases for MODIS

Total absorption optical depth by atmospheric gases for SGLI

0.06

0.06

channel ch06
channel ch07

0.05

0.04

τabs

τabs

0.04

0.03

0.03

0.02

0.02

0.01

0.01

0

channel SW03
channel SW04

0.05

0

2

4

6

8
Atmospheric level (km)

10

12

14

2

4

6

8
Atmospheric level (km)

10

12

14

Figure 7: Gaseous absorption profile comparison between two SWIR channels. Left: MODIS
channel 6 and 7; Right: SGLI channel SW03 and SW04

and one “non-absorptive” channel from SGLI for the cloud masking over snow. As shown in
Fig. 9, ice is non-absorbing in SGLI channel SW01 and “slightly” absorbing in SGLI channel
SW03 and SW04. So we followed our strategy to select two pairs of SWIR channels as
SW01+SW03 and SW01+SW04 and did simulation accordingly. The results are consistent
with our expectation as shown in Fig. 11, which demonstrates that cloudy and cloud-freecases are distinguishable. We actually used a similar channel choice in Chen et al. [2014],
where MODIS channels 5 and 7 were used with Aqua MODIS data, and its performance wasvalidated by collocated CALIOP measurements.

SGLI Algorithm Description

15

2.1 C1: Cloud mask and surface classification

1

1
SGLI
SW01

MODIS
Ch05

MODIS SGLI
Ch06
SW03

SGLI
SW01

MODIS SGLI
Ch07 SW04

MODIS
Ch05

MODIS SGLI
Ch06
SW03

MODIS SGLI
Ch07 SW04

0.95

0.85

Asymmetry parameter

Single scattering albedo

0.9← r =100µm
eff

0.8

0.75

0.7

r =400µm →

0.95

eff

← reff=100µm

0.9

0.65
reff=400µm →

0.6

0.55
0.85

0.5

1

1.2

1.4

1.6

1.8

Wavelength (µ m)

2

2.2

2.4

1

1.2

1.4

1.6

1.8

2

2.2

2.4

Wavelength (µ m)

Figure 8: Inherent optical properties of snow particles (grain size 100 µm and 400 µm) in
the SWIR region. Left: Single-scattering albedo; Right: Asymmetry parameter.

−3

x 10
2

1.4
MODIS
Ch05

MODIS SGLI
Ch06 SW03

MODIS SGLI
Ch07 SW04

1.35

1.5

1.3

1

1.25

0.5

1.2

1

1.2

1.4

1.6

1.8

Wavelength (µ m)2

2.2

Imaginary Part

Real Part

SGLI
SW01

0
2.4

Figure 9: Refractive index of ice in the SWIR region. Solid line: real part; Broken line:
imaginary part.

SGLI Algorithm Description

162.1 C1: Cloud mask and surface classification

1

Single scattering albedo (MODIS Ch7/SGLI SW04)

snow
ice cloud, solid column
ice cloud, hollow column
ice cloud, aggregate
ice cloud, rosette−4
ice cloud, rosette−6
ice cloud, plateice cloud, dendrite

reff = 10µm
reff = 15µm
reff = 20µm
reff = 10µmreff

reff = 40µm

= 15µm0.95
reff = 20µm

reff = 60µm

reff = 40µm
0.9

r

eff

= 100µm

reff = 60µm

0.85
0.85

0.9

0.95

1

Single scattering albedo (MODIS Ch6/SGLI SW03)

Figure 10: Single-scattering albedo of snow and ice cloud particles in the two SWIR channels.
Results for MODIS channels 6 and 7 are shown in green and for SGLI channels SW03 andSW04 in blue.

SGLI Algorithm Description

172.1 C1: Cloud mask and surface classification

0.35

0.35

snow surface with different grain sizes
pure ice clouds (τi = 0.1,0.5,1.0,2.0,5.0)

0.25

0.2

0.15

0.1

0.05

0
0.1

snow surface with different grain sizes
pure ice clouds (τi = 0.1,0.5,1.0,2.0,5.0)

0.3

SGLI channel SW4 (2.20 µm) reflectance

SGLI channel SW3 (1.63 µm) reflectance

0.3

0.25

0.2

0.15

0.1

0.05

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

SGLI channel SW1 (1.05 µm) reflectance

0
0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

SGLI channel SW1 (1.05 µm) reflectance

Figure 11: The same as Figure 6 but for alternative SGLI channels. Left: SGLI SW01 andSW03; Right: SGLI SW01 and SW04.

SGLI Algorithm Description

18

2.2 C2: snow, sea ice and atmospheric parameter retrievals

2.2

C2: snow, sea ice and atmospheric parameter retrievals

Basically three key steps are used to retrieve the snow and atmospheric parameters, (i) aforward radiative transfer model, (ii) a neural network training algorithm and (iii) a nonlinear
optimal estimation method. Figure 12 shows a flow chart of the retrieval algorithm for snow
and atmospheric parameters.

SGLI Algorithm Description

19

2.2 C2: snow, sea ice and atmospheric parameter retrievals

2.2.1

Forward radiative transfer model

The discrete-ordinate radiative transfer (DISORT)
model for the atmosphere-snow system was used
to calculate the radiance at the top of the atmosphere as a function of atmospheric and snow
physical parameters [Stamnes et al., 1988b, 2007b].
This model has a pseudo-spherical treatment forsolar beam attenuation in a curved spherical-shellatmosphere [Spurr, 2002], and it has been validated against Monte-Carlo results [Gjerstad et al.,2003]. We calculate the extinction coefficient,
single-scattering albedo, and phase function usinga Mie code based on conventional light-scattering
theory, under the assumption of spherical grain
shapes.

START
Linearized
DISORT
Radiance data
! (!, st, rs; ")NN training

Satellite
Data: ! (")#

Inverse model

Updatep(input)Snow model In the snow model, since we know
the size distribution of snow particles [Aoki et al.,2000], an effective snow grain radius is considered
for the optical properties of snow. The size distribution of snow particles is a log-normal size distribution with a geometric standard deviation of 1.6
measured by Grenfell & Warren [1999a] in Antarctica. The refractive indices of ice are taken fromthe data compiled by Warren & Brandt [2008] (Fig.
13). For a mixture of snow/ice and impurities, theeffective optical properties of the mixture are obtained by weighting each component so that thevolume extinction coefficient β̃ext , the single scattering albedo ω and the phase function asymmetry
factor g of an external mixture are obtained as follows:!"!
!"
Converge?
#$%!
!, st
st, rsr(output)

END

Figure 12: Flow chart of the atmospheric and snow parameter retrieval
algorithm. All notations are describedin the text.

imp
imp
ice
β̃ext = (1 − f )β̃ext+ f β̃ext= [(1 − f )ρs κice
ext + f ρimp κext ],

ω =
g =

imp imp(1 − f )ω ice κice
κext
ext + f ω,
imp
ice
(1 − f )kext + f kextice
imp
(1 − f )ω ice κice
+ f ω imp κimp
ext g
ext g
,
imp
(1 − f )κice
+
f
κ
sca
sca

(2)
(3)
(4)

where κext and κsca are the mass extinction and scattering cross sections, respectively, andthe superscripts ice and imp denote the components for ice and impurity, respectively, f is
the mass-fraction of the impurity in unit [ppmw], ρs is the snow mass density, and ρimp is
the impurity mass density. For snow impurities, we used a black carbon model [Hess et al.,SGLI Algorithm Description20

2.2 C2: snow, sea ice and atmospheric parameter retrievals1998] for the refractive index of the snow impurity (Fig. 13). The snow optical depth foreach homogeneous layer within the snowpack isice,` `ice `τs` = β̃exth = ρ`s βext
h /ρice

and

ice,`
βext
=

`
)
3Q̃`ext (reff
,
`
4reff

(5)

`
, ρ`s and h` are the effective snow grain radius, the snow mass density and thewhere reff`
) is the extinction
thickness of the `th layer. ρice is the density of pure ice, and Q̃`e (reffefficiency.

Figure 13: Real and imaginary part of the refractive index of ice [Warren & Brandt, 2008]and black carbon [Hess et al., 1998].
In the previous snow algorithm, we assumed that the snow grains have spherical shapes.
A comparison of BRDF measurements obtained by the Cloud Absorption Radiometer (CAR)
flown on a NASA aircraft over Elson Lagoon in Barrow, Alaska, see Fig. 14, shows that using
non-spherical (Voronoi) particle shapes to describe snow grains in our coupled atmospheresurface radiative transfer model (CRTM) provides a better match with the BRDF measurements than using spherical particles. For the best match, the grain size and impurity are
50 µm and 0.001 PPMW, and the aerosol optical depth is 0.01. These values are consistent with those derived from MODIS data closest in time. This BRDF comparison provides
an opportunity to test which scattering phase function yields the best agreement with themeasured BRDFs. We found that use of the phase function computed from an assemblyof non-spherical Voronoi particles yields a better match with the measured BRDFs than aHenyey-Greenstein phase function with asymmetry factor obtained from Mie computations.
We have found that the BRDF produced by an assembly of spherical particles cannot match
SGLI Algorithm Description

21

2.2 C2: snow, sea ice and atmospheric parameter retrievals

the measured one, no matter how we adjust the grain size, the impurity concentration or theaerosol properties. This predicament implies that the HG phase function is unsuitable forsimulations of the BRDF of snow. Therefore, in year 2015, we updated our snow model touse Voronoi IOPs in order to improve the performance of the snow retrieval algorithms forrealistic snow surfaces. The Voronoi particle IOPs has been provided by Dr. Aoki’s groupusing the ellipse ratio M03. Comparisons of retrieved snow grain sizes from MODS dataobtained using spherical and Voronoi particle models will be discussed in Section 3.2.2.
Phase function for Voronoi particles In many RTMs an expansion of the phase function in Legendre polynomials is used. This approach requires accurate computation of the
expansion coefficients (phase function moments), which is challenging for phase functions
with sharp forward peaks occurring for scattering by particles that are large compared to
the wavelength. Experimenting with phase functions produced by an assembly of nonspherical Voronoi particles (Ishimoto et al., 2010), we found that special techniques areneeded to properly treat such a phase function. Due to the very strong peak in the forward
scattering direction, it is very difficult to get the correct phase function moments, neededin the radiative transfer calculations, by numerical integration, and, if handled improperly
significant errors will occur. In order to deal with this problem, we experimented with thedelta-fit method of Hu et al. (2000) and extended it to treat the entire phase matrix. As an
example, we show in Fig. 15 all phase matrix elements for a collection of spherical particles
computed using a Mie code, and approximated by the delta-fit method. The truncation of
the forward peak by the delta-fit method is quite pronounced. Figure 16 shows I, Q, U and
V components of the reflected light at the TOA computed by our vector radiative transfer
code (VDISORT) using the exact and the delta-fit approximated phase matrices. The good
match for all output angles indicates that our delta-fit approximation to the phase matrix
elements yields reliable TOA Stokes parameters for reflected light. We then applied thedelta-fit method to the phase matrix generated from an assembly of Voronoi particles. As
shown in Fig. 17 the phase matrix resulting from application of the delta-fit method closely
matches the original phase matrix. Hence, our delta-fit treatment of the Voronoi phasematrix appears to be robust.
Surface elevation correction In cryosphere areas, many places are in high elevation regions, such as Greenland and Antarctica. In these areas, the near surface atmosphere gas
distribution will be different from that at sea level. The gas absorption and molecular scattering layer will be thiner, when the elevation is higher. In order to correct the elevation
effect, we simulated satellite radiances for different elevations (0–4 km). Then using a neuralnetwork technique, we developed a height correction table to convert the satellite measured
radiance from any elevation (0–4 km) to the corresponding sea surface radiance. By doing
so the snow/ice parameter retrieval algorithm, based on the sea level radiances, could be applied for the retrieval. Figures 18 and 19 show the difference between uncorrected elevation
radiance and the sea level radiance. The difference is up to 10-15% at short wavelengths andup to 40% in the NIR, due to the strong Rayleigh scattering at short wavelengths and thestrong water vapor absorption in the NIR. The smallest difference is in the 0.865 µm and1.050 µm channels, since there is weaker absorption and less scattering. So the elevation
SGLI Algorithm Description

222.2 C2: snow, sea ice and atmospheric parameter retrievalscorrection is very important in our snow retrieval algorithm. Figure 20 shows a comparison
between corrected sea level radiance and the true sea level radiance. Our correction performance is pretty good. The error in the corrected radiance is within about 0.3%.

Snow layer depth Our algorithm was designed with 2 snow layers to explore the snow
vertical structure. Instead of the critical depth used previously, we have fixed the depth of
the first snow layer to be 1.5 cm and the second layer to be 98.5 cm. This new snow first layer
depth is in the range of 1.64 µm channel penetration depth. Since the 1.64 µm channel has
strong snow absorption, and a low S/N ratio, we changed our algorithm to use the 1.05 µm
channel instead of the 1.64 µm channel for snow parameter retrieval. An additional benefit
is that the 1.05 µm channel is less affected by humidity and surface elevation change than
the 1.64 µm channel. In section 3.2, we will apply the new algorithm to MODIS data.
Sea ice model The inherent optical properties (IOPs) of snow can be parametrized interms of the snow grain size and impurity concentration. Similarly, the IOPs of ice dependon scattering/absorbing inclusions (primarily air bubbles and brine pockets) and impurities
embedded in the ice [Jin et al., 1994; Hamre et al., 2004; Jiang et al., 2005; Stamnes et al.,2011]. If we can retrieve information about the size of the snow grains and scattering ice
inclusions as well as impurity concentrations, we may use a radiative transfer model for the
coupled atmosphere-snow-ice-ocean system (hereafter referred to as CRTM) to compute thesnow/ice BRDF and albedo. Assuming that snow grains and ice inclusions can be adequately
represented by a collection of spherical particles, we may write the IOPs, i.e. the absorption
and scattering coefficients and the scattering phase function as
αp (λ) =Z rmaxrmin

2πr Qα (r)n(r)dr ,R rmaxpp (λ, Θ) =rmin

σp (λ) =Z rmaxrmin

πr2 Qσ (r)n(r)dr ,pp (λ, Θ, r)n(r)dr,rmin n(r)dr

R rmax(6)(7)

where the absorption or scattering “efficiency” Qα (r) or Qσ (r) is defined as the ratio of theabsorption or scattering cross section for a spherical particle of radius r to the geometricalcross section πr2 , and n(r) is the particle size distribution. For a specific value of theparticle radius r, we can compute Qα (r), Qσ (r), and pp (λ, Θ, r) using Lorenz-Mie theory,
but the integrations in Eqs. (6)-(7) require knowledge of the particle size distribution n(r)which is usually unknown. Equations (6)-(7) can be considerably simplified by making thefollowing assumptions [Hamre et al., 2004; Stamnes et al., 2011]: (i) The particle distribution
is characterized by an effective radius
R rmaxn(r)r3 dr
2rmin n(r)r dr

ref f = Rrrmin
maxwhich obviates the need for an integration over r. (ii) The particles are weakly absorbing,so that
16π ref f mi,p 1[m3rel − (m2rel − 1)3/2 ](8)Qα (r) ≈3λmrel
SGLI Algorithm Description

232.2 C2: snow, sea ice and atmospheric parameter retrievals

where mi,p is the imaginary part of the refractive index of the particle, λ is the wavelength
in vacuum, and mrel = mr,p /mr,med is the ratio of the real part of the refractive index of
the particle (mr,p ) to that of the surrounding medium (mr,med ). (iii) The particles are large
compared to the wavelength (2πr/λ >> 1) which implies Qσ (r) = 2. The scattering phasefunction may be represented by the one-parameter Henyey-Greenstein (HG) phase function,
which depends only on the asymmetry factorg ≡ hcos Θi =1Z 1p(Θ) cos Θ d(cos Θ).
2 −1With these assumptions, Eqs. (6)-(7) become:1[1 − (m2rel − 1)3/2 ]fV ,mrel
1 − g2pp (λ, Θ) =.
(1 + g 2 − 2g cos Θ)3/2
αp (λ) = α(λ)σp (λ) =3 fV,
2 ref f

(9)(10)

Here α(λ)
is the absorption coefficient of the material of which the particle is composed, and
4π R
3fV ≡ 3 n(r)r3 dr ≈ 43 πreff ne , where ne = number of particles per unit volume with radius
ref f . As described by Stamnes et al. [2011] this approach can be extended to work very well
for wavelengths less than about 2.8 µm.
Aerosol model For aerosol retrieval, we assume a fixed aerosol model in the SGLI snow
retrieval algorithm as we did for GLI. But we do consider real humidity information in
the aerosol model. Different geographical regions are expected to have different aerosol
properties, and all aerosol models will be obtained from the OPAC compilation [Hess et al.,1998]. So only the aerosol optical depth will be retrieved. For the Arctic and Antarctica, wewill directly use the OPAC Arctic and Antarctic aerosol models. For mid-latitude regions,
there are three types of aerosol models available in OPAC (Continental Average, Continental
Clear, Continental Polluted). We are not going to retrieve the relative humidity (RH),
instead we will directly use a given (near-real time) RH value (obtained from other satellitedata) in our aerosol retrieval. We pre-calculated the IOPs for each aerosol type with RH= 0%, 50%, 70%, 80%, 90%, 95%, 98%, and 99%. Then we linked these 8 humidity values
to provide a continuum of aerosol IOPs by using spline interpolation on humidity. The
humidity is considered to be a given variable which is used in the forward radiative transfer
calculation as well as in our retrieval algorithm. Figure 21 shows the impact of humidity
on the TOA radiance. We compare the TOA radiance for RH = 0% humidity with the RH
values from 50% to 99%. From Fig. 21, we note that in the visible the RH effect could be 4
- 11%, while in the NIR, this effect is smaller, and about 1 - 4%. These results demonstratethat use of real time RH information will be very important for the snow retrieval becauseof the large impact on radiances, especially in the visible spectral range.
In our current algorithm, the aerosol model will be selected based on the geographicalarea [Hess et al., 1998]. The near-real time humidity will be obtained from another sensor,
and used in the retrieval code as a given parameter. Thus, the snow and aerosol retrieval
will be based on near-real time humidity information. The atmospheric model will be the
US standard model.
SGLI Algorithm Description242.2 C2: snow, sea ice and atmospheric parameter retrievals

2.2.2

Neural Network (NN) Training

In order to retrieve the snow and atmospheric parameters simultaneously, we need to use an
inversion method. The traditional lookup table (LUT) interpolation method is not suitable
for simultaneous and accurate retrieval of multiple parameters. Instead, we use a nonlinear, iterative optimal estimation method [Li et al., 2008]. But this approach is very time
consuming because we need to call the forward model repeatedly to compute radiances andJacobians required in the nonlinear, iterative optimal estimation method described in thenext section. In order to solve this problem, we replaced the forward radiative transfer modelwith a radial basis function neural network (RBF-NN) to establish a relationship between
the retrieval parameters (referred to as P ) and the spectral radiances (referred to as R). A
two-layer RBF-NN was employed here for training of the interpolating radial basis functions.
(a) Spectral radiance to parameters training (R2P training): We start by traininga RBF-NN in which the spectral radiances and the solar zenith angles are used as input data,
whereas the snow and aerosol parameters are the output data. The benefit of this training isthat we can use the results directly to estimate the initial snow and atmospheric parameters
for the retrieval. Thus, by using the training coefficients, we can directly calculate the snow
and atmospheric parameters given by
Pk =

N
X

akj exp{−b2

j=1

Nr
X

(Ri − cji )2 } + dk ,

(11)

i=1

where N is the total number of neurons, and Nr is the sum of the number of input radiances
and the number of solar zenith angles. The Ri ’s are the input parameters, and akj , b, cji anddk are the optimized coefficients.
(b) Parameters to spectral radiance training (P 2R training): The purpose of thisRBF-NN is to provide a relationship between the retrieval parameters and the radiances anduse this analytic relationship to derive an analytic expression for the Jacobians. Both theradiances and the Jacobians are required in the nonlinear optimal estimation described in
the next section. The input data are the snow and atmospheric parameters (i.e. snow grainsize, snow impurity concentration, and aerosol optical depth) and the solar zenith angles.
The output data are the radiances in all SGLI channels. The spectral radiance Ri in channel
i obtained from the radial basis neural network (RBF-NN) is given as:
Ri =

N
X
j=1

aij exp{−b2

NpX(Pk − cjk )2 } + di ,

(12)

k=1

where N and Np are the total number of neurons and input retrieval parameters, respectively.
The Pk ’s are the input parameters, and aij , b, cjk and di are the optimized coefficients. Notethat the radial basis functions functions in Eqs. (11) and (12) are trained on the logarithmsof Ri and Pk , and that normalization has been used to adjust the magnitude of the different
input parameters, because the magnitudes of grain size and optical depth are quite different.
SGLI Algorithm Description25

2.2 C2: snow, sea ice and atmospheric parameter retrievalsThis difference could decrease the accuracy of the Gaussian radial basis functions [Eq. (12)].
The normalization alleviates this problem.
The Jacobians (partial derivatives) are also required for the retrieval by nonlinear optimal
estimation as discussed in the next section. These Jacobians can be calculated by analyticaldifferentiation as follows
Jik = −2NX

aij b2 (Pk − cjk ) exp{−b2

j=1

2.2.3

Np
X

(Pl − cjl )2 }.

(13)

l=1

Nonlinear optimal estimation

We use cost function minimization techniques appropriate to nonlinear iterative spectral
fitting [Li et al., 2008]. The retrieval parameters are assembled into a state vector x =[P1 , P2 , . . . , Pn ]T , and the measured radiances into another vector ymeas = [R1 , R2 , . . . , Rm ]T .
The update of the state vector xi at iteration step i is given by [Rodgers, 2000]:
xi+1 = xa + Gi ymeas − yi + Ki (xi − xa )

h

i

(14)

Gi = Sa KTi (Ki Sa KTi + S )−1 .

(15)

The measurement vector ymeas has covariance matrix S corresponding to noise , yi =F(xi , b) are simulated radiances generated by the forward model F(xi , b) which is a (nonlinear) function of xi and b (vector of model parameters not retrieved but sources of error).
Ki is the Jacobian matrix of simulated radiance partial derivatives with respect to xi . The apriori state vector is xa , with covariance Sa . Gi is the gain matrix of contribution functions.
The inverse process starts from an initial guess x0 ; often set to xa . A convergence criterion
checks progress towards the solution that minimizes the cost function. Besides measurement
noise and a priori covariance (smoothing error), other sources of error are uncertainties inthe elements of b, and forward model uncertainties due to physical and/or mathematical
simplifications.
In order to speed up the retrieval, our algorithm used the RBF-NN method instead of
the traditional LUT approach. This new algorithm has several benefits:• We use the information available in all channels to solve the inverse problem by the
optimal estimation method to produce simultaneous retrieval of all desired parameters.
• We use the P 2R neural network to compute the simulated spectral radiance RN N andJacobian J N N , which speeds up the retrieval by a factor of 100 or more compared to
the LUT appraoch.
• We use the R2P neural network to obtain a good first guess of our retrieval. Thisapproach will effectively decrease the number of iterations required for the optimal
estimation algorithm to converge.

SGLI Algorithm Description26

2.2 C2: snow, sea ice and atmospheric parameter retrievals

2.2.4

Snow and sea ice surface albedo retrieval

In the current algorithm, we provide 2 types of broadband snow albedos and one type ofsea-ice broadband albedo in 3 wavelengths ranges (VIS: 0.3 − 0.7 µm; NIR: 0.7 − 2.8 µm;
SW: 0.3 − 2.8 µm). For snow covered surfaces, we have 2 ways to compute the broadband
albedo: direct albedo and indirect albedo. For sea-ice covered areas, we provide only the
direct sea-ice albedo.
Snow indirect albedo From the physical parameters inferred from our snow/ice retrieval
algorithms and the corresponding IOPs derived from them, we may use the CRTM to calculate the broadband albedo, the surface BRDF and the spectral albedo of the coupled system.
These snow/ice parameters include the snow/ice particle size and the impurity concentration. The broadband snow/ice albedo obtained in this manner is called indirect albedo. It
cloud be the black (transparent) sky (if we ignore the atmosphere) or the actual or blue-sky
snow/ice albedo (if we include the atmosphere and the retrieved aerosol optical depth), but
it is calculated without invoking the Lambertian assumption because the full BRDF of thesurface is accounted for in the CRTM. This snow/ice albedo estimate, based on retrieved
snow/ice IOPs, provides results that are more accurate than those produced by the standard
MODIS broadband albedo algorithm [Li et al., 2007]. In this algorithm, we provide only
the broadband albedo for SGLI standard products, but we have the potential to provide the
BRDF and spectral albedo too.
Snow direct albedo Although the indirect albedo provide a lot of important information,
one key limitation of this indirect albedo estimation is that its accuracy will depend on theaccuracy of the retrieved snow parameters. If the retrieved snow grain size and impurity
concentrations are inaccurate, the inferred albedo results will be inaccurate too. So how do
we validate the retrieved snow grain size and impurity concentration? One option would beto rely on field measurements, but one is limited by the sparsity of field data, which may beavailable only for a few days of the year at a few locations. Another limitation is that fielddata do not allow for validation of snow retrievals on a pixel by pixel basis. In view of thisspace-time sparsity of field measurements and the critical need for validations, we provide
another albedo, which called the direct albedo.
The direct albedo approach consists of retrieving the actual snow surface albedo directlyfrom the measured satellite radiances. Using our CRTM for the coupled atmosphere-snow
surface system, we can simultaneously simulate the satellite radiances and the snow surface
albedo as a function of snow grain size and impurity concentration and aerosol optical depth.
These simulations can be used in conjunction with a neural network to estimate directly thesurface albedo and the aerosol optical depth from 8 SGLI channels (or 7 MODIS channels).
In this manner one obtains the actual snow surface albedo of the coupled atmosphere-snow
surface system. Because this actual snow surface albedo is inferred directly from the measured satellite radiances, it is not influenced by inaccuracies in our retrieved snow parameters. This new direct snow albedo algorithm, can be used for fast and accurate broadband
snow albedo estimation, and the difference between the old (indirect) and the new (direct)
albedo estimations will provide a valuable indicator about the quality of the retrieved snow
SGLI Algorithm Description27

2.3 C2: Surface temperature retrieval

parameters.
The quality of the retrieved snow parameters depends to a large extent on the performance of the cloud mask and the surface classification algorithm. Inferior quality of theretrieved snow parameters in a given pixel is mainly due to imperfect identification of thinclouds, aerosols, melting snow, ice (instead of snow), or partial vegetation (forest) cover, andlack of removal of these effects in the TOA radiances. For such “fake” snow pixels, which arefrequently encountered in satellite images, the TOA radiances do not depend on just puresnow, and the performance of the snow parameter retrieval algorithms is mostly determinedby such “fake” snow pixels. To alleviate this problem, we are providing both direct and indirect albedo estimates. By comparing these two albedo estimates we will be able to provide
a snow retrieval quality flag for every pixel in an image. Figure 22 is an example of how we
can use the difference between the two albedo values to estimate the quality of the retrieved
snow parameters. We note that there are many pixels around the edge of Greenland withlarge relative differences, implying that the retrieved snow parameters are untrustworthy.
These areas mostly correspond to the sea-ice or melting snow.
Sea-ice direct albedo In year 2014, we used the CRTM to model the radiative transfer
in the coupled atmosphere-sea ice system [Stamnes et al., 2011]. The satellite measured
radiance over sea-ice was simulated as a function of the abundance of air bubbles and brine
pockets in the sea ice as well as its thickness. The direct albedo was computed in a similar
way as described for the snow direct albedo. A neural network technique was used to link thesimulated satellite radiance and pre-calculated values of three broadband surface albedos.
The simulated radiances were computed for sea ice conditions ranging from new ice, first year
ice to multiple year ice. Table 3 lists the sea ice types and ranges used in the simulations.
Each of the three sea ice types includes 10,000 simulation cases, which gives total of 30,000
cases for the neural network training. The direct sea ice albedo model (synthetic) data yields
a very good retrieval as can be seen in Fig. 23. In next year we will test the sea-ice albedo
retrieval algorithm using MODIS data.
Table 3: Sea-ice types and ranges used in the simulations
item
sea-ice thickness (m)air bubble volume fractionbrine pocket volume fraction

2.3

New iceFirst-year ice Multiple-year ice
0.05 - 0.3
0.3 - 1.0
1.0 - 3.0
0.16 - 0.35
0.11 - 0.24
0.01 - 0.15
0.005 - 0.012 0.008 - 0.04
0.012 - 0.07

C2: Surface temperature retrieval

Based on the GLI snow surface temperature retrieval algorithm, we developed a similar
algorithm for the MODIS and the SGLI sensors. A split-window technique was adopted toestimate the surface temperature in the polar regions by using MODIS channels 31 (11 µm)
and 32 (12 µm) as well as SGLI channels T1 (10.8 µm) and T2 (12 µm). This algorithm
SGLI Algorithm Description28

2.3 C2: Surface temperature retrieval

consists of two parts: one is using model emissivities to generate the algorithm, another isbased on field-measurements [Hori et al., 2006] to generate the algorithm. Both the sensor’s
viewing angle and the snow grain size were taken into account in this algorithm. Thisalgorithm could be applied not only to snow-covered sea and land surfaces, but also to a
mixture of snow/ice and melt-ponds. It works only under clear-sky conditions. Even though
the technique used to create this algorithm is similar to that used in estimating sea and land
surface temperatures in general, this algorithm has been developed specifically for the polar
regions and for use with SGLI /MODIS measurements.
2.3.1

Background

The Arctic is particularly sensitive to global climate change. Accurate estimate of surface
temperature could provide an early signal of climate change. The surface temperature in thepolar regions controls sea-ice growth, snow melt, and surface-atmosphere energy exchange.
During the past decade, significant progress has been made in estimation of sea surface temperature and snow/ice surface temperature from satellite thermal infrared data. A common
approach for estimating surface temperature is to relate satellite data to surface temperature observations with a regression model. A radiative transfer model can be applied tomodel satellite radiances and brightness temperature. A large set of atmospheric profiles
and measured surface temperature will be used to simulate the satellite measurements under a wide range of conditions. The success of the algorithm depends primarily upon thevariability of surface and atmospheric characteristics. In order to correct for atmosphericeffects, split-window techniques centered around 11 and 12 µm are commonly employed.
Such approaches have been used for SST and IST retrieval [Minnett, 1990; Llewellyn-Joneset al., 1984; Barton, 1985; Key & Haefliger, 1992; Key et al., 1997; Wan & Dozier, 1996;
Wan, 2008] and they have been applied to AVHRR and ATSR (Along Track Scanning Radiometer) data using two “split-window” infrared channels at approximately 11 and 12 µm.
For the SST algorithm an absolute accuracy of 0.5-1 K (root-mean-square error or RMSE)
has been obtained [Llewellyn-Jones et al., 1984; McClain et al., 1985]. IST accuracies of
0.3-2.1 K relative to measured or modeled surface temperatures have been reported by Key
et al. [1997] using the split-window technique. Land surface temperature (LST) estimationis generally less accurate due to the larger variability of surface conditions [Price, 1983; Wan& Dozier, 1996].
The surface emissivity is defined as the ratio of the actual radiance emitted by a given
surface to that emitted by a black body at the same kinetic temperature. It is a key
parameter in this split-window technique. The thermal emissivity of snow/ice has an angularand spectral dependency [Dozier & Warren, 1982; Warren, 1982; Hori et al., 2006; Hall et al.,2008]. Also it varies with snow grain size especially at larger grain sizes [Hori et al., 2006;
Hall et al., 2008]. The most significant improvement in this version, except for the effectof the sensor angle, is that we made the snow emissivity depend on grain size. The fieldemissivity data were collected in Japan and Alaska during the year 2002-2004 [Hori et al.,2006]. We keep both the model emissivity version and the field-emissivity version here forresearch purposes. The model version includes only the sensor angle effect. The coefficients
listed in Table 9 are based on the MODIS channels. Since response functions for the SGLISGLI Algorithm Description29

2.3 C2: Surface temperature retrieval

channels are now available, in year 2015, we updated the coefficients based on the response
functions for the SGLI channels (see Table 4 - Table 8).
2.3.2

Algorithm description

In this study the MODTRAN radiative transfer model is employed to simulate the radiances
measured by the satellite sensor (e.g., SGLI) using the directional snow emissivities thatare computed with the DISORT radiative transfer model. Since MODTRAN has a 2 cm−1
spectral resolution, it is accurate enough for the purpose of this study. To simulate radiancesin SGLI and MODIS thermal channels, daily temperature and humidity profiles are usedin the MODTRAN radiative transfer model. Radiosonde ascents over the entire Arctic aretaken from the NCEP/NCAR (National Center for Atmospheric Research) Arctic Marine
Rawinsonde archive. This data set contains 17,659 reports of ship (marine) rawinsondes (i.e.,radiosondes: tracked from the ground by radar to measure variations in wind direction andwind speed with altitude) for the region north of 65◦ . Its record extends from 1976 to 1996.
Sounding data from this NCAR Rawinsonde archive cover different atmospheric conditions(such as those caused by regional and seasonal variations).
The simple approach for atmospheric correction is to measure radiation from a given fieldof view at two or more window frequencies having different atmospheric absorption. The
surface temperature can then be estimated as a linear combination of measured brightness
temperatures at these frequencies:
TS = a(θ) +nXbi (θ)Ti ,

(16)

i=1

where TS is the measured surface temperature, θ is the satellite scan angle, a(θ) and bi (θ) arescan angle-dependent coefficients, Ti is the measured brightness temperature in each thermalchannel i, and n is the total number of channels used. The minimization of errors in the TS
measurements relies on correct choice of the coefficients a(θ) and bi (θ).
There are two ways to determine the coefficients. One approach is to relate satellite
observations to surface temperature measurements with a simple regression model. However,for a robust solution a relatively large set of high-quality in-situ temperature and satellite
data is required. The other approach is the simulation method. A radiative transfer model isused together with a large set of atmospheric profiles to simulate the satellite measurements
under a wide range of atmospheric conditions and surface temperatures. The simulated
measurements are then used with a set of assigned surface temperature values to derive thecoefficients, again by regression analysis.
Instead of computing a different set of coefficients for each scan angle increment, as shownin (16), here we use the equation
TS = a + bT11 + c(T11 − T12 ) + d[(T11 − T12 )(sec θ − 1)],(17)

where TS is the estimated surface temperature (in K), T11 and T12 are the brightness temperatures (in K) at 11 µm (SGLI channel T1, MODIS channel 31) and 12 µm (SGLI channel T2,
MODIS channel 32), respectively, and θ is the sensor scan angle. The coefficients a, b, c, and
SGLI Algorithm Description30

2.3 C2: Surface temperature retrievald are derived from multilinear regression. This approach using two “split-window” infrared
channels at approximately 11 µm and 12 µm is commonly employed for surface temperature
retrieval [Minnett, 1990; Llewellyn-Jones et al., 1984; Barton, 1985; Key & Haefliger, 1992;
Key et al., 1997; Wan & Dozier, 1996].
To ensure high accuracy the MODTRAN radiative transfer model is employed for simulating the radiances measured by SGLI/MODIS. In MODTRAN a narrow-band model isused for computing gaseous optical depth from the HITRAN database with wavenumber
steps of 1 cm−1 , which implies a nominal spectral resolution of 2 cm−1 at FWHM (FullWidth at Half-Maximum). Multiple scattering is also included in the radiative transfer
model by combining MODTRAN and DISORT. For the SGLI and MODIS instruments,
since the sensor scan angle lies between 0◦ and 50◦ for SGLI and between 0◦ and 65◦ for
MODIS, our simulations are done for viewing angles in the range of 0◦ − 65◦ . The built-in
standard subarctic winter and summer atmospheric profiles including trace gases and thebackground aerosol model in MODTRAN are used in our simulations. Blanchet & List [1983]
showed that the volume extinction coefficient of Arctic haze is generally of the same orderof magnitude as for troposphere aerosols. Thus, we use the troposphere background aerosol
model instead of Arctic haze. The calculations of the retrieval coefficients in Eq. (17) using
both SGLI and MODIS channels are presented in this document for testing and validating
the algorithm. The sensor response functions both for SGLI channels T1 and T2 and for
MODIS channels 31 and 32 are used to compute radiances at the top of the atmosphere, andthe radiances are then converted to brightness temperature by use of the Planck function.

Table 4: The SGLI coefficients in Eq. (17) for model emissivity algorithm.
Temperature range (K)≤ 240
240 -260
260-270
270 - 275
≥ 275a-0.9420168-1.700981-0.5846105-3.221689
2.843076b
1.003281
1.006895
1.003292
1.012690
0.9904238

c2.080047
1.668042
1.329147
1.455035
1.562278d0.2917113
0.4842514
0.5522773
0.4839154
0.4033772

Corr. coef.
0.9991578
0.9998550
0.9975470
0.9915375
0.9963182

Table 5: The SGLI coefficients in Eq. (17) for fine dendrite snow in field emissivity algorithm.
Temperature range (K)< 240
240 -260
260-270
270 - 275

SGLI Algorithm Descriptiona-1.090729-1.788184-0.5331097-3.630045

b1.004445
1.007761
1.003598
1.01469931c2.084182
1.656945
1.317247
1.441280d0.3006721
0.5195864
0.5939672
0.5123382

Corr. coef.
0.9991575
0.9998645
0.9975812
0.9904801

2.3 C2: Surface temperature retrieval

Table 6: The SGLI coefficients in Eq. (17) for medium granular snow in field emissivity
algorithm.
Temperature range (K)< 240
240 -260
260-270
270 - 275a-1.248099-2.110133-0.7415222-4.578181

b1.005447
1.009418
1.004917
1.018801c
2.083595
1.751293
1.331195
1.415663d0.2596135
0.4588805
0.6152086
0.5270104

Corr. coef.
0.9991561
0.9998040
0.9980019
0.9942210

Table 7: The SGLI coefficients in Eq. (17) for coarse grain snow in field emissivity algorithm.
Temperature range (K)< 240
240 -260
260-270
270 - 275a-1.299243-2.174245-0.7264020-4.302524b1.005852
1.009710
1.004918
1.017884c
1.927811
1.758209
1.348086
1.413851d0.2369569
0.4572233
0.6092463
0.5349259

Corr. coef.
0.9991557
0.9998040
0.9980387
0.9943376

Table 8: The SGLI coefficients in Eq. (17) for sun crust snow in field emissivity algorithm.
Temperature range (K)< 240
240 -260
260-270
270 - 275a-0.7152216-1.733492-1.223238-4.154361b1.000973
1.005832
1.004883
1.015466c2.055423
1.758956
1.373944
1.419622d0.1783278
0.3113386
0.4460161
0.4283974

Corr. coef.
0.9991505
0.9998055
0.9981163
0.9959655

Table 9: The MODIS coefficients in Eq. (17) for model emissivity algorithm.
Temperature range (K)≤ 240
240 -260
260-270
270 - 275
≥ 275

SGLI Algorithm Descriptiona-1.624761-2.019964-5.224606-2.013436-0.4194403

b1.008296
1.009724
1.022082
1.009982
1.00408732c
2.800785
2.500067
1.568301
1.558308
1.821280d-0.9120480-1.009879
0.1110692-1.298285
1.644374

Corr. coef.
0.9991578
0.9998550
0.9975470
0.9915375
0.9963182

2.3 C2: Surface temperature retrieval

Figure 14: BRDF comparison between snow model simulation (voronio and spherical particle) and CAR measurement over Elson Lagoon, Barrow on April 7, 2008. The modeled
BRDF was computed at height 635m, which is the aircraft height. Top 3 panels are the3600 BRDF. lower 4 panels are BRDF at principal plane, 450 , 1350 plane and cross plane,
separately.
Table 10: The MODIS coefficients in Eq. (17) for fine dendrite snow in field emissivity
algorithm.
Temperature range (K)< 240
240 -260
260-270
270 - 275

SGLI Algorithm Descriptiona-1.793135-2.072019-4.873211-1.887228

b1.009592
1.010481
1.021244
1.01003833c2.802395
2.503243
1.713263
1.624779d-0.8154156-0.9555640-0.3119795
0.9791106

Corr. coef.
0.9991575
0.9998645
0.9975812
0.9904801

2.3 C2: Surface temperature retrieval

2o

Mie scattering matrix elements, original vs delta−fitted, Kokhanovsky (2010) aerosol, rg = 0.3µ m, [ln(sigmag)] = = 0.8464, λ = 0.412µ m, ΘT = 5 , NMOM = 150

F11(phase function)

4

F33/F1110

1.5
originaldelta−fitoriginaldelta−fit

3

10

1

2

10

0.5

1

10

0

0

10

−0.5

−1

10

−1

−2

10

0

45

90

−1.5

135

0

45

F12/F1190

135

F34/F110.2

0.2
originaldelta−fitoriginaldelta−fit

0.1

0.1
0
0

−0.1
−0.2

−0.1

−0.3
−0.2

−0.4
−0.5

−0.3

−0.6
−0.4
−0.7
−0.5

0

45

90

−0.8

135

0

45

90

135

Figure 15: Elements of the Stokes scattering matrix for an ensemble of spherical aerosol
particles: Original versus delta-fit approximation.

SGLI Algorithm Description34

2.3 C2: Surface temperature retrieval

φ=0°

φ=90°

0.4

0.05

0.3

0.04

0.2

0.03

φ=180°
0.1
0.08
0.06
0.04

0.1
0
0

0.02

20

40

60

80

0.02

0.01
0

φ=0°

20

−3

0

6

−0.01

x 10

40

60

80

Mie
Delta fit

0
0

20

φ=90°

40

60

80

60

80

60

80

60

80

φ=180°
0.01

5

0

4

−0.01

3

−0.02

−0.02
Mie
Delta fit

−0.03
−0.04
−0.05
0

20

40

60

80

2
0

20

°

40

60

80

−0.03
0

°

φ=0

φ=90

1

−8

0

4

−0.002

0.5

20

x 10

40
°

φ=180

Mie
Delta fit

2

−0.00400−0.006

−0.5
−1
0

−2

−0.008
20

40

60

80

−0.01
0

φ=0°

20

−6

1

3

x 10

40

60

80

φ=90°

2

x 10

−4

−1
−0.5

60

80

−3
0

Mie
Delta fit

−6

−2
40

φ=180°

−2

0

20

40

0

1

−1
0

20

−12

2

0.5
0

−4
0

20

40

60

80

−8
0

20

40

Figure 16: Stokes parameters of reflected light at TOA computed using exact Mie and deltafit approximated phase matrices. From top to bottom, I, Q, U , and V components. The
exact Mie results were obtained using 300 moments whereas 150 moments were used toobtain the delta-fit results.

SGLI Algorithm Description

35

2.3 C2: Surface temperature retrieval

Voronoi(snow) scattering matrix elements, original vs delta−fitted, reff = 50µ m, λ = 0.38µ m, ΘT = 1o, NMOM = 200

F11(phase function)

6

F22/F1110

F33/F1118
originaldelta−fit5

105
originaldelta−fit

16

0

14

4

originaldelta−fit

10

12

−5

3

10

10

2

10

8

1

6

10

−10

−15

4

0

10

2

−20

−1

10

0

−2

10

0

50

100

150

−2

0

50

F44/F11

100

150

0

50

F12/F11

6

100

150

F34/F111

4

originaldelta−fit4−25

originaldelta−fit0.5

originaldelta−fit
3

2
0

0

2

−2

−0.5

−4

1
−1

−6

0

−8

−1.5

−10

−1

−2
−12
−14

0

50

100

150

−2.5

0

50

100

150

−2

0

50

100

150

Figure 17: Elements of the Stokes scattering matrix for an ensemble of non-spherical Voronoi
particles: Original versus delta-fit approximation.

SGLI Algorithm Description36

2.3 C2: Surface temperature retrieval

Figure 18: The simulated satellite radiance relative difference (%) between different elevations and sea level.
380 nm443 nm50

50
RMSE: 1.2263Eï02

RMSE: 5.6653Eï03

40

40

40

40

35
30
25
20
15
10

0
ï10

35
30
25
20
15
10
5

R

ï5

TOA

0

5

relative difference (%)

0
ï10

10

Case No. percentage (%)

45

Case No. percentage (%)

45

Case No. percentage (%)

35
30
25
20
15
10
5

R

ï5

TOA

0

5

relative difference (%)

865 nm

0
ï10

10

R

TOA

30
25
20
15
10

5

0
ï10

10

40

20
15
10
5
0
ï10

Case No. percentage (%)

45

40

Case No. percentage (%)

45

40

25

25
20
15
10
5

ï5

0

5

RTOA relative difference (%)

10

0
ï10

35
30
25
20
15
10
5

ï5

0

5

RTOA relative difference (%)

10

0
ï10

10

5

10

2210nm45

30

5

RMSE: 1.2344Eï01

40

35

0

relative difference (%)

50
RMSE: 3.9728Eï02

45

30

ï5

TOA

50
RMSE: 2.1657Eï03

35

R

1640nm50
RMSE: 2.4983Eï03

0relative difference (%)

1050nm

50

35

5

ï5

Case No. percentage (%)

Case No. percentage (%)

670 nm50
RMSE: 2.0285Eï0245

5Case No. percentage (%)

530 nm50
RMSE: 2.9252Eï02

45

35
30
25
20
15
10
5

ï5

0

5

RTOA relative difference (%)

10

0
ï10

ï5

0

RTOA relative difference (%)

Figure 19: The histogram of simulated satellite radiance difference between uncorrected
elevation radiance and sea level radiance.
SGLI Algorithm Description37

2.3 C2: Surface temperature retrieval

380 nm443 nm

10

5

20

15

10

5

0

RTOA relative difference (%)

865 nm15

100

RTOA relative difference (%)

0
ï0.5

0.5

0

RTOA relative difference (%)

10

5

20

15

10

5

0

RTOA relative difference (%)

0.5

0
ï0.5

2210nm

25

20

15

10

5

0

RTOA relative difference (%)

0.5

0.5

RMSE: 7.7766Eï04

25

Case No. percentage (%)

Case No. percentage (%)

15

0

RTOA relative difference (%)

30
RMSE: 5.0045Eï04

25

20

10

0
ï0.5

0.5

30
RMSE: 1.4696Eï0425

151640nm30
RMSE: 1.6477Eï0420

51050nm

30

0
ï0.5

20

5

0
ï0.5

0.5

25

Case No. percentage (%)

0ï0.5

RMSE: 2.7154Eï04

25

Case No. percentage (%)

Case No. percentage (%)

Case No. percentage (%)

15

30
RMSE: 5.0664Eï0425

20

670 nm30
RMSE: 8.0266Eï04

25

Case No. percentage (%)

530 nm30
RMSE: 1.1010Eï03

Case No. percentage (%)

30

0
ï0.5

20

15

10

5

0

RTOA relative difference (%)

0.5

0
ï0.5

0

RTOA relative difference (%)

0.5

Figure 20: The histogram of simulated satellite radiance difference between corrected sea
level radiance vs. the true sea level radiance.

Relative difference from RH=0% on TOA radiance

12
RH=50%
RH=70%
RH=80%
RH=90%
RH=95%
RH=98%
RH=99%

10

8

6

4

2

0
400

600

800

1000
1200
Wavelength(nm)

1400

1600

1800

Figure 21: Humidity effect (%) on TOA radiance. Aerosol model: Arctic aerosol from OPAC.
Solar zenith angle 50◦ , sensor zenith angle 65◦ , relative azimuth 120◦ , aerosol optical depth
0.08 at 859 nm, first layer grain size 80 µm, second layer grain size 800 µm, impurity 0.01
ppmw.
SGLI Algorithm Description

38

2.3 C2: Surface temperature retrieval

Figure 22: The relative albedo differences between the Terra and Aqua images obtained onJune 27, 2012 over Greenland.

Figure 23: Sea-ice albedo model data retrieval.
SGLI Algorithm Description

392.3 C2: Surface temperature retrieval

Figure 24: Flow chart of the surface temperature retrieval algorithm.
For the model emissivity version, the directional surface emissivities of snow are calculated based on an extended version of the DISORT radiative transfer model appropriate forthe coupled atmosphere-surface system. This algorithm allows us to study the bidirectionalreflectance and directional emissivity for a surface covered with different types of snow andsea ice. The coefficients a, b, c, and d in Eq. (17) for SST retrieval are calculated within thefollowing four temperature ranges: (i) T11 ≤ 240 K; (ii) 240 K ≤ T11 ≤ 260 K; (iii) 260 K≤ T11 ≤ 270 K; and (iv) 270 K ≤ T11 ≤ 275 K.
Table 11: The MODIS coefficients in Eq. (17) for medium granular snow in field emissivity
algorithm.
Temperature range (K)< 240
240 -260
260-270
270 - 275a-0.9379274-2.202930-5.629201-2.846899

b1.004896
1.010195
1.023626
1.012921c
2.737998
2.145490
1.206063
1.494790d-1.335196-0.4138538-1.213855
1.672925

Corr. coef.
0.9991561
0.9998040
0.9980019
0.9942210

For the field emissivity version, the field measured emissivity ([Hori et al., 2006]) wasused to calculate the coefficients a, b, c, and d. There are 4 different types of snow according
to the snow grain size. They are the fine dendrite snow; medium granular snow, coarse grainsnow, and sun crust. When the grain size increases, the emissivity decreases. We used the
SGLI Algorithm Description40

2.3 C2: Surface temperature retrieval

Table 12: The MODIS coefficients in Eq. (17) for coarse grain snow in field emissivity
algorithm.
Temperature range (K)< 240
240 -260
260-270
270 - 275a-1.206548-2.377483-5.616658-2.792477b1.006264
1.011157
1.023869
1.013001c
2.743953
2.087973
1.192855
1.489832d-1.425086-0.2680590
1.248157
1.701098

Corr. coef.
0.9991557
0.9998040
0.9980387
0.9943376

Table 13: The MODIS coefficients in Eq. (17) for sun crust snow in field emissivity algorithm.
Temperature range (K)< 240
240 -260
260-270
270 - 275a-1.338066-2.513868-5.590907-3.277824b
1.007215
1.012142
1.024266
1.015255c
2.738751
2.012193
1.125157
1.484559d-1.453996-0.1153396
1.502457
1.794687

Corr. coef.
0.9991505
0.9998055
0.9981163
0.9959655

same four temperature ranges in the model emissivity version for these four snow grain types.
So in this version, we first establish to which snow type the retrieved first layer grain size
belongs. The second step is to check the brightness temperature at 11 µm, to decide which
temperature group to use. The last step is to use Eq. (17) to calculate the temperature.
For temperatures larger than 275 K, we assume that a mixture of snow/ice and meltponds occurs. Because we have no field measurement data for this situation, we will use
model emissivity. In such cases the surface emissivities are assumed to be an area-weighted
sum of snow and water emissivities. The weights are 0.2 for snow and 0.8 for water. Obviously, this approximation could lead to an uncertainty in the estimated surface temperature.
The procedure used to determine the coefficients in Eq. (17) may be summarized asfollows (see Fig. 24). The surface physical properties, observed atmospheric profiles, sensor
response functions, and snow directional emissivity data, calculated using the DISORT RT
model or taken directly from field measurements, are input into MODTRAN to simulate
radiances at SGLI/MODIS channels for a wide range of atmospheric conditions. In thesesimulations, the temperature of the first layer just above the surface is taken to be the surface
temperature. The simulated radiances are integrated with the sensor response functions for
the SGLI/MODIS channels:
Ri =

Z λ2
λ1

R(λ)φi (λ)dλ,(18)where Ri is the simulated radiance, and φi is the sensor response function for channel i.
Then, the integrated, simulated radiances Ri are converted to brightness temperatures Ti .
A least-squares multilinear regression is used to determine the coefficients in Eq. (17). The
coefficients in Eq. (17) for the MODIS channels are presented in Tables 9-13, respectively.
The correlation coefficients between estimated and actual surface temperature are also given
SGLI Algorithm Description

41

2.4 Aerosol over snow retrieval by use of polarization channels

in Tables 9-13.
In order to develop an algorithm using the split-window technique specifically for the polar regions, atmospheric profiles from the NCEP/NCAR Arctic Marine Rawinsonde Archive
are used for simulating the sensor-measured radiances. More than 7,000 radiosonde profiles are used for the surface temperature algorithm development. The data set from the
NCEP/NCAR Arctic Marine Rawinsonde Archive covers a large range of atmospheric conditions including different seasons and different locations across the Arctic.

2.4

Aerosol over snow retrieval by use of polarization channels

Aerosols have a major impact on the radiative energy balance and climate. In view ofthe scarcity of cryospheric field measurements, satellite remote sensing has a great advantage
by providing accurate estimation with good temporal and wide spatial coverage. However,owing to the high brightness of snow surface that dominates the TOA radiance, and thelarge variation in snow reflectance determined by grain size (snow age and temperature) andembedded impurities, aerosol retrieval over high albedo snow surfaces remains a challenging
problem.
In year 2015, we attempted to address this problem by using the information in SGLI
non-polarized channels available for snow retrieval to obtain a direct simultaneous fitting of
the measured TOA radiances to derive an optimum set of aerosol and snow properties. For
this purpose we used a linearized radiative-transfer model for the coupled atmosphere-snow
system, in combination with a non-linear, iterative optimal estimation (OE) inverse method
for simultaneous retrieval of aerosol and snow properties. In the validation work, we foundthat the aerosol retrieval result is not as good as desired, because it is difficult to separate thecontribution of aerosols to the TOA radiance from that due to snow impurities. Therefore,
we tried to add SGLI polarized channels in our simultaneous retrieval to improve the aerosol
results, since the polarized signal will be more sensitive to the small aerosol particles than tothe large snow particles. We already did some preliminary work using our Vector Radiative
Transfer Model (VDISORT: ?; ?) to compute the polarized radiance I = [I, Q, U, V ]t at the
TOA. For this purpose we used a standard aerosol model [Hess et al., 1998] and the IOPs of
the snow grains were computed either from Mie theory (assuming spherical particles) or byusing an assembly of non-spherical Voronoi particles [?]. The snow impurities were assumed
to consist of soot particles internally mixed with the snow grains.
Figure 25 shows the sensitivity of TOA scalar reflectance
RI = πI/µ0 F0and polarized reflectanceq

RP = π Q2 + U 2 /µ0 F0in the SGLI channel at 0.6735 µm to variations in the soot concentration for a snowpackconsisting of Voronoi particles with a fixed grain size of 100 µm. The aerosol model used
was OPAC summertime (Hess et al., 1998).
Figure 26 shows a similar sensitivity test for variations in aerosol optical thickness. As
expected, the results show that for large snow grains the polarized reflectance is insensitive
SGLI Algorithm Description

422.4 Aerosol over snow retrieval by use of polarization channels∆φ = 0◦∆φ = 90◦∆φ = 180◦

0.6

0.55

Reflectance

0.65
0.6

0.5

0.5

0.55

0.45

0.5

0.4
0.4

0.45
0.4

0.3

0.35

0.35

Polarized Reflectance

0

20

40

60

80

0.2
0

0.11

0.16

0.1

0.14

0.09

0.12

0.08

0.1

0.07

0.08

0.06

0.06

20

40

60

80

0

20

40

60

80

0.08
soot=0.01
soot=0.1
soot=1.0
soot=10.0

0.06
0.04
0.02
0.05
0

20

40

60

80

0.04
0

20

40
60
View Angle (°)

80

0
0

20

40

60

80

Figure 25: Reflectances for several values of the soot concentration in the snow. Upper
panels: scalar reflectance RI . Lower panels: Polarized reflectance RP . The aerosol model
was from OPAC summertime (Hess et al., 1998) with an optical thickness of 0.05. The solar
zenith angle was θ0 = 45◦ .
to the soot concentration, whereas the polarized reflectance in the forward direction (∆φ ∈[0◦ , 90◦ ]) is sensitive to aerosol optical thickness. From this sensitivity study, we can conclude
as follows:• It is not possible to distinguish between AOT (aerosol optical thickness) effects andabsorption by snow impurities using radiance-only observations.
• The AOT can be distinguished from snow impurity by using the polarized reflectance.
• The polarized reflectance in the forward direction (∆φ ∈ [0◦ , 90◦ ]) is more sensitive to
AOT.
The SGLI polarized sensor is designed to be tilted into the forward scatting direction (”Forward Looking” in the NH and ”Backward Looking” in the SH) (see Fig. 27). Thus, it provides
a good opportunity to improve aerosol retrieval over bright snow surfaces by using SGLI’s
polarization channels. We will continue this work in the next 3 years. We will also explore
the advantage of combining the scalar radiance measurements obtained from the other SGLI
channels appropriate for retrieval of aerosol/snow properties with the polarized radiances
obtained from the two SGLI channels at 0.6735 µm and 0.8685 µm. This approach will also
SGLI Algorithm Description43∆φ = 0◦∆φ = 90◦

Reflectance

0.9
0.8

0.8

0.75

0.7

0.7

0.6

0.65

0.5

0.6

0.4

0.55

0.7
τa = 0.05
τa = 0.10
τa = 0.20

0.6
0.5
0

20

40

60

80

0.16
Polarized Reflectance∆φ = 0◦

0.14

0

20

40

60

80

0.5
0

0.25

0.2

0.2

0.15

0.15

0.1

0.1

0.05

20

40

60

80

20

40

60

80

0.12
0.1
0.08
0.06
0

20

40

60

80

0.05
0

20
40
60
View Angle (°)

80

0
0

Figure 26: Reflectances for several values of the aerosol optical thickness. Upper panels:
scalar reflectance RI . Lower panels: Polarized reflectance RP . The aerosol model was from
OPAC summertime (Hess et al., 1998 and the snowpack consisted of an ensemble of 100 µm
Voronoi particles externally?? mixed with 1.0 ppm black carbon. The solar zenith angle
was θ0 = 45◦ .
allow us to use information obtained from multiple angles, which is expected to not only
improve the quality of retrieved snow properties, but also help solve the problem associatedwith the aerosol retrievals over bright surfaces.

3

Algorithm Application and Validation

3.1

Application and validation of the C1 algorithm

We tested our methodology on MODIS data and compared our results with corresponding
results produced by the latest version (Collection 6) of the MODIS Cloud Mask (MOD35 and
MYD35 for MODIS Terra and Aqua, respectively). In addition to the image-based tests we
did a direct comparison of C1 and MYD35 with the CALIOP 1 km cloud layer product for awhole year of data collected over Greenland. This long-term statistical validation using theactive cloud detection scheme of CALIOP provides the most reliable assessment currently
available of both the C1 and the MODIS cloud masks, and more information can be obtained

SGLI Algorithm Description

44

3.1 Application and validation of the C1 algorithm

Figure 27: The SGLI polarization sensor orbit direction.

SGLI Algorithm Description

45

3.1 Application and validation of the C1 algorithm

by comparing the results from different seasons. Data from year 2007 were chosen since thisyear is the first year for which CALIOP has a whole year of observations over Greenland,since CALIOP onboard the CALIPSO satellite was launched April 28, 2006. Only land
pixels were included in the statistics since we want to focus on the Greenland land area. The
blue bands of MODIS Terra are continually degrading at a faster rate than those of MODIS
Aqua. We therefore want to use old data not influenced by this degradation to avoid its
impact on C1’s NDSI snow identifier. From June 1 to August 10 we used the Sub-Arctic
Summer atmospheric constituent profiles [Anderson et al., 1986] in the dynamic thresholdcalculation whereas the Sub-Arctic Winter profiles were chosen for the rest of the year. In
the tests we configured C1 to use different SWIR channels for the MODIS sensors on Terra
and Aqua:
1. MODIS channel 6 (1.64 µm) and channel 7 (2.13 µm) were used in the test with MODIS
Terra data. Since these two channels are close to SGLI channels SW3 (1.63 µm) andSW4 (2.20 µm) we consider tests with MODIS Terra data to be a simulation of C1’s
performance on the future SGLI sensor.
2. MODIS channel 5 (1.24 µm) and channel 7 (2.13 µm) were used with MODIS Aqua
data due to the detector problem of MODIS Aqua channel 6. MODIS Aqua data can be
collocated with CALIOP measurements. Hence, we can use CALIOP as a benchmark
to evaluate the performance of our algorithm. Also it would be interesting to see how
it performs using different SWIR channels and we consider it a good way to test our
methodology.
3.1.1

Image-based tests over Greenland

As a cloud mask designed to work in the polar regions, one of the main targets of theC1 algorithm is that it should deliver reliable cloud mask results over Greenland, which
is one of the most difficult places for a cloud mask algorithm because of the high altitudeof the Greenland ice sheet and the extremely low surface temperature. Figure 28 shows acomparison between C1 and MOD35 over Greenland with the northwest coast of Greenland(Qaanaaq and Peary land) located at the upper left corner of each panel. In the color
composite, the 0.65 µm band is assigned to the red, the 1.6 µm band to the green and the12.0 µm band to the blue. The clouds are usually in yellow (low-level water clouds), white(mid-level water clouds), orange or red (thick ice clouds with low cloud top temperature)and sometimes blue (thin ice clouds). Snow areas are usually shown in magenta. In thecloud mask results (bottom panels) clouds with cloudy confidence levels from low (darkgray, 66% < clear confidence level ≤ 99%), medium (gray, 33% < clear confidence level ≤66% ) to high (white, clear confidence level < 33%, mostly cloudy) are plotted together
with snow (red), sea ice (purple), open land (green) and open ocean (blue). For MODIS
Cloud Mask results the “confidently cloudy” and “probably cloudy” categories are treatedas cloudy scenes and “confidently clear” and “probably clear” are treated as clear scenes. In
the circled part of the image areas with thick ice clouds are apparent from either the RGB
or the 11 µm BT plots (very low temperatures). MOD35 (lower left) mis-classified this area
as snow, while C1 (lower right) correctly masked out this cloudy area. Our tests indicate
SGLI Algorithm Description

46

3.1 Application and validation of the C1 algorithm

that in the polar regions the MODIS cloud mask will sometimes miss some of the cirrus
clouds or multi-layer clouds with very low temperature. The drastic temperature change
during daytime and temperature inversions may create problems for the BT and BTD tests,
whereas the reflectance-based scheme of our algorithm is largely unaffected by temperatureeffects.
3.1.2

Direct comparison with CALIOP using MODIS Aqua data

The launch of NASA’s Cloud Aerosol Lidar and Infrared Pathfinder Satellite Observation(CALIPSO) satellite with its onboard CALIOP instrument provides vertically resolved cloud
and aerosol information. When accurately collocated with MODIS Aqua measurements
CALIOP can provide valuable evaluation of the cloud mask performance of C1 and MYD35
from an active cloud detection perspective. Liu et al. [2010] used CALIOP data to evaluate
the accuracy of MODIS Cloud Mask products over Arctic sea ice and found that the MODIS
Cloud Mask generally performs better over open ocean than over Arctic sea ice surfaces. In
our study we are focusing on the Arctic region and on the Greenland plateau, in particular.
The collocated data from CALIOP and MODIS Aqua was taken from March to Octoberof 2007 and the collocation method presented by Holz et al. [2008] was used. We chose tostart our comparison from March and end by October because C1 is a daytime algorithm
and there is very little data in the Winter season. We are taking the results from CALIOP
as the ‘truth’ to assess the performance of C1 and MYD35. Some detailed statistics of the
hit rate and the Hanssen-Kuipers Skill Score or True Skill Score (TSS) [Hanssen & Kuipers,
1965] from three time periods in 2007 are shown in Table 14. The hit rate is defined as:
Hit rate =Ncld,hit + Nclr,hit
Ntotal(19)

where Ncld,hit is the number of pixels for which C1 or MYD35 agreed with CALIOP thata cloud was detected, and Nclr,hit is the number of pixels for which they agreed that clearsky was detected and Ntotal = Ncld,hit + Ncld,miss + Nclr,hit + Nclr,miss is the total number of
MODIS Aqua pixels collocated with CALIOP measurements, including pixels for which thetwo algorithms disagreed. The TSS defined as:
TSS =(Ncld,hit · Nclr,hit − Ncld,miss · Nclr,miss )(Ncld,hit + Ncld,miss ) · (Nclr,hit + Nclr,miss )(20)

is widely used to evaluate the effectiveness of a prediction. It provides useful informationof the hit rate with respect to the false alarm rate (miss rate), and will remain positive
as long as the hit rate is higher than the false alarm rate. A higher TSS means a morereliable prediction and identification from the algorithm. Figure 29 shows a comparison of
hit rates and TSS obtained for C1 and MYD35 using CALIOP 1 km cloud layer results as abenchmark. From the statistics we see that C1’s performance is generally better than that of
MYD35 especially from late-Spring to mid-Autumn with higher hit rate and TSS score. In
early Spring (March, April)/late Autumn (October) C1 and MYD35 perform very similarly,
but C1 generally has a higher hit rate but slightly lower TSS, indicating a higher false alarm
rate. We checked images from the above seasons and found that our algorithm would produce
SGLI Algorithm Description

473.1 Application and validation of the C1 algorithm

280

270

260

250

240

230

220

40

40

35

35

600

600 30

30

800

800

1000

1000

200

200

400

400

25

25

20

20

15

15

1200

1200

1400

1400

10

10

1600

1600

1800

1800 5

5

2000

2000
200

Clouds'

400

600

800

Snow'over'land'

1000

200

1200

Sea'ice'Open'water'

400

600

800

1000

1200

Open'land'

Figure 28: C1 vs. MOD35 Collection 6 (MOD35 C6) for a Terra MODIS image, obtained on
July 29, 2005 over Greenland. Upper left: False color RGB composite using MODIS band 1(0.65 µm) for red, band 6 (1.6 µm) for green, and band 32 (12 µm) for blue. Upper right: 11µm brightness temperature, note the thick ice clouds with very low temperature (≈ 220-240K) in circled areas. Lower left: The MOD35 C6 cloud mask. Lower right: The C1 cloudmask.

SGLI Algorithm Description

48

3.1 Application and validation of the C1 algorithm

more false clouds compared to the middle part of the year. Since the atmosphere over theGreenland plateau can become very dry with large variations in specific humidity during
Spring and late Autumn [Ettema et al., 2010], use of a fixed atmospheric profile (Sub-Arctic
Winter) may account for the higher false cloudy rate in these seasons. However, a small
grain size due to surface frost or diamond dust near the snow surface [Grenfell & Warren,
1999b], and the effect of surface roughness (sastrugi) under large solar zenith angles [Warren
et al., 1998; Kuchiki et al., 2011], may also contribute to the high false cloudy rate.
100
MYD35 Hit Rate
C1 Hit Rate
MYD35 Hanssen Kuiper
C1 Hanssen Kuiper

85

80

95
90
85
80

Hit Rate (%)

75
75

70
65

70

60
55

65

50

Hanssen Kuiper Skill Score (%)

90

45
60

40
35

55

2

3

4

5

6

7

8

9

10

30
11

Month

Figure 29: C1 and MYD35 Collection 6 (MYD35 C6) vs. CALIOP over Greenland, 2007.

Table 14: C1 and MYD35 vs CALIOP in Greenland, 2007
Time period

Total
MYD35
matchups HRa(%)
Mar. 01 to May 15 535283
69.34
May 15 to Aug. 23 1291045 74.66
Aug. 23 to Oct. 31 289405
73.69
Overall
2115733 73.18

3.1.3

C1
HR
(%)
70.64
76.90
74.34
74.96

MYD35
TSSb
(%)
40.71
49.97
42.21
47.08

C1
TSS
(%)
42.03
53.78
42.44
49.76

Cross comparison of C1 using MODIS Terra and Aqua data

As mentioned above, the C1 algorithm is designed for SGLI implying that tests using MODIS
Terra data would be the best way to assess its performance on data to be collected with the
SGLI sensor since both sensors have the 1.64 µm band available. However MODIS Terra
cannot be collocated with CALIOP so there is no direct way to evaluate the performanceof C1 using the 1.64 µm channel. As a substitute we compared the cloud fraction overGreenland derived from C1 using both MODIS Aqua and Terra data. If the cloud fractions
SGLI Algorithm Description493.1 Application and validation of the C1 algorithm

derived by C1 from the two sensors are close at the same time period of the year, then we
may consider it to be an indirect validation of the performance of C1 on MODIS Terra data.
Figure 30 shows the 4-day averaged (using data from day 1-4, 5-8, 9-12, etc.) cloud fraction
over Greenland in 2007. The cloud fraction is calculated as: fcld = Ncld /(Ncld + Nclr ) where
Ncld is the total number of pixels classified as cloudy and Nclr is the number of pixels classified
as clear by the algorithm during that period. It can be seen that the cloud fractions derived
by C1 from the two sensors are generally very close throughout the entire year, indicatingthat the performances of C1 using the 1.64 µm channel on MODIS Terra and the 1.24 µm on
MODIS Aqua are comparable. Larger difference is observed in the Spring season as Terra’s
cloud fraction is lower than that of Aqua by about 5%, which may be related to the highfalse cloudy rate obtained when using MODIS Aqua data as discussed in Section 3.1.2. Since
the higher degree polynomials needed for interpolation of the 1.24 µm channel on Aqua may
introduce larger error in our algorithm, we consider the results obtained from Terra data to
be closer to the actual cloud fraction in the Spring season.
1
C1 on MODIS Terra data
0.95

C1 on MODIS Aqua data

0.9

Cloud fraction

0.85

0.8

0.75

0.7

0.65

0.6

0.55
50

100

150

200

250

300

Day of year

Figure 30: C1 derived 4-day averaged cloud fraction over Greenland, 2007.

3.1.4

Summary

Multi-spectral satellite instruments like AVHRR, MODIS and VIIRS provide high spatial
resolution observations over the entire Arctic region. Their data products can be used toderive useful information such as cloud fraction and distribution, snow/ice coverage, surface
temperature etc. which will help us better understand Arctic climate and its variability. An
accurate cloud mask over the Arctic region is critical in order to obtain reliable products
from satellite data. The traditional methods include thermal IR based cloud tests (BT6.7 ,BTD3.7−11 etc.) and reflectance-based cloud tests (employing combinations of VIS/SWIR
tests with simple solar/viewing geometry correction). Water vapor bands like the 7.2 µm
band and the carbon dioxide band at 14.2 µm are also applied to provide additional information. In some of these studies radiative transfer calculations were performed to determine
SGLI Algorithm Description50

3.1 Application and validation of the C1 algorithm

the cloud screening thresholds by simulating different kinds of clouds for a variety of different
atmospheric and surface conditions in the Arctic.
The C1 algorithm is developed specifically for the SGLI sensor, which has very few
thermal IR channels available. A cloud screening scheme that mainly uses the reflectanceof two SWIR channels is established. Compared to other reflectance-based cloud screening
methods the thresholds of our method are dynamically determined by the solar/viewing
geometry for each satellite pixel based on comprehensive radiative transfer calculations thattake into account (snow) surface BRDF, surface elevation and the atmospheric profiles ofscattering/absorbing molecules. Compared to the traditional thermal IR based methods ourscheme has the following advantages:
1. The thresholds are less sensitive to the drastic temperature changes and frequent temperature inversions occurring in the Arctic.
2. The reflectance thresholds are based on computed snow surface reflectance rather thancloud properties since it is difficult to simulate every possible cloud configuration that
may occur under complex atmospheric and surface conditions in the Arctic.
3. A smaller number of satellite channels is required for desirable performance. Sensors
with no water vapor and carbon dioxide channels should be able to benefit from acloud mask similar to C1.
The validation of our cloud masking scheme was performed using data from MODIS
Terra, Aqua and CALIOP over Greenland in 2007. The image-based tests show that ouralgorithm improves the detection of cold ice clouds and cloud edges compared to the thermal
IR based MODIS Cloud Mask algorithm. Statistical validation with CALIOP on MODIS
Aqua data shows that our C1 cloud mask generally has a higher hit rate and a higher TSS
score than the MODIS cloud mask. The cloud fraction derived using MODIS Terra and Aqua
data shows that the two configurations of our algorithm using different SWIR channels (1.64µm and 2.13 µm for MODIS Terra data, 1.24 µm and 2.13 µm for MODIS Aqua data) leadto consistent performance. More work is needed for the early Spring and late Autumn overGreenland since a higher false cloudy rate is observed. Appropriate profile information of
the atmospheric absorption for each SWIR channel may be required because a fixed model
atmosphere is expected to provide an inadequate representation of the atmospheric conditions over Greenland in all seasons. This problem may be addressed by using surface specifichumidity data from the AWS sites on Greenland to correct the standard model atmosphere.
Results from the water vapor channels of MODIS would also be helpful. Another possible source of error might be the presence of optically thin ice clouds (τ = 0.1 - 0.5) overGreenland in these seasons which presents a challenge to C1 and other reflectance-based
cloud mask algorithms. As indicated by our model calculations the reflectances of theseclouds depend strongly on the ‘background’ surface reflectances and detailed simulations using well-defined surface properties (like snow with different grain sizes) might be a good way
to proceed. Other possible sources of mis-classification include surface frost, diamond dust,
sastrugi and surface slope. Data from the VIIRS instrument, which has SWIR channels very
similar to those of SGLI, and which can be collocated with CALIOP observations, would be
SGLI Algorithm Description

51

3.2 Application and validation of the C2 algorithm

very useful for further assessment of the performance of our C1 algorithm. Our existing test
and validation framework can easily be adjusted for application to VIIRS data.

3.2

Application and validation of the C2 algorithm

3.2.1

Comparison with GLI algorithms and field measurements

We compared the SGLI algorithm results with the GLI algorithm for the same MODIS
images. All these MODIS images corresponded to field measurement data. Aoki et al.
[2007a] describe in some the detail the snow pit work for all the field data. For each snow
pit work site, we calculated all the satellite pixels within 1 km circle, and then we used
the averaged retrieval data of these pixels to compare with the field measurements. The
purposes of our work here are (1) to validate our new SGLI algorithm including retrieval
feasibility and accuracy, (2) to compare with GLI results to highlight improvements, and (3)to discuss how to improve the next version of this algorithm. Figure 31 shows the snow grain
size comparison. The SGLI algorithm retrieval results are shown in the top panels and GLI
results in the bottom panels. Both first layer (left panels) and second layer (right panels)
snow grain sizes were compared with the field measured data. The field measurements havethe measured grain size range, so we plotted the measurement range here for comparison.
From Fig. 31 we see that SGLI has more data, which means that the SGLI algorithm retrieved
all available MODIS pixels around the field measurement site. But for the GLI algorithm
the retrieval failed for some of these pixels. Thus, the SGLI algorithm has a much better
retrieval capability than the GLI algorithm. This superior performance could be due to thefact that the SGLI algorithm is based on a nonlinear inverse method which uses the radiancein all channels simultaneously. We used different colors to indicate the data for different snow
pit sites, as well as wet snow and dry snow marked by “*” and “o”. Most retrieved snow
grain sizes lie in the measurement range, but the performance is better for dry snow thanfor wet snow. For wet snow, actually that is the case in Nakashibetsu. As we discussed in
Section 4.2, due to high reflectance from a “sun crust” surface, we got underestimated grain
size for the top layer and overestimated grain size for the bottom layer.

SGLI Algorithm Description

52

3

10

Barrow
Nakashibetsu
SaromaAbashira
BihoroRetrieved Rs2 (µ m) ï SGLIRetrieved Rs1 (µ m) ï SGLI

3.2 Application and validation of the C2 algorithm

2

10

1

10
1
10

2

10

Measured R

10
1
10

3

10

at 0ï0.5 cm (µ m)Retrieved Rs2 (µ m) ï GLI

Retrieved Rs1 (µ m) ï GLI

10*: wet snowo: dry snow

1

2

10

s1

3

10

2

10

1

10
1
10

3

2

10

3

10

Measured Rs2 at 0.5ï5 cm (µ m)

3

10

2

10

1

2

10

10
1
10

3

10

Measured Rs1 at 0ï0.5 cm (µ m)

2

10

3

10

Measured Rs2 at 0.5ï5 cm (µ m)

Figure 31: Snow grain sizes retrieved from the SGLI algorithm (top) and the GLI algorithm(bottom) compared with field measurements. Left: first layer grain size versus 0 − 0.5
cm layer field measurements. Right: Second layer grain size versus 0.5 − 5 cm layer fieldmeasurements.

SGLI Algorithm Description

53

0

10

ï1

10

ï2

10

ï2

10

ï1

10

0

10

1

10

1

10

0

10

ï1

10

ï2

10

ï2

10

ï1

10

0

10

1

10

Measured total carbon (ppmw)*: wet snowo: dry snow

0

10

ï1

10

ï2

10

ï2

0

10

10

1

10

0

10

ï1

10

ï2

10

ï2

10

0

10

Measured total snow impurity (ppmw) Measured total snow impurity (ppmw)

Retrieved impurity Cs (ppmw) ï GLI

Retrieved impurity Cs (ppmw) ï GLI

Measured total carbon (ppmw)

1

10

Retrieved impurity (ppmw) ï SGLI (dust)

Barrow
Nakashibetsu
SaromaAbashira
Bihoro

1

10

0

10

ï1

10

ï2

10

ï2

0

10

10

Retrieved impurity (ppmw) ï GLI(dust)

1

10

Retrieved impurity Cs (ppmw) ï SGLIRetrieved impurity Cs (ppmw) ï SGLI

3.2 Application and validation of the C2 algorithm

1

10

0

10

ï1

10

ï2

10

ï2

10

0

10

Measured total snow impurity (ppmw) Measured total snow impurity (ppmw)

Figure 32: Snow impurity concentrations retrieved from SGLI algorithm (top) and GLI
algorithm (bottom) compared with field measurements. Left: retrieved impurity versus
measured total carbon. Middle: retrieved impurity versus measured total impurity. Right:
equivalent dust impurity vs. measured total impurity.
Figure 32 shows the snow impurity comparison. As for the snow grain size, the SGLI
results are shown in the top panels, and the GLI results in the bottom panels. In the snow
model, the snow impurity concentration is assumed to be the same in both snow layers. We
compare with the field measured averaged (0 - 5 cm) impurity data, and the snow impurity
is assumed to consist exclusively of black carbon. In the in-situ snow pit work, the main
composition of snow impurity was found to be mineral dust [Aoki et al., 2007a]. Since
the spectral dependence of these two components are different in the visible, the resultscannot be compared directly. In Fig. 32, left panels, we compare the retrieved impurity
with the measured total black carbon part only, and we note that there is an overestimationof the retrieved impurity concentration. If we compare with the total measured impurity
concentration which include both black carbon and mineral dust (middle panels of Fig. 32),
the retrieved impurity shows an underestimation. Considering the difference between mineral
dust and black carbon, and noting that mineral dust has about 10 times the effect of black
carbon, we may convert the retrieved soot (black carbon) amount to an equivalent dust
amount (right panels of Fig. 32), and then we get good agreement with the total measured
impurity concentration. The SGLI and GLI algorithms yield similar results except that SGLI
algorithm gets more retrievable pixels. We also noticed that for wet snow over Nakashibetsu,
SGLI Algorithm Description

543.2 Application and validation of the C2 algorithm

Figure 33: Snow parameter retrievals: Terra MODIS images over Greenland, Day 170, 2005.
the retrieved snow impurity is questionable too. We will test different snow impurity types
and mixing conditions in our next version of the impurity retrieval algorithm, as well as thewet snow condition.
3.2.2

Application to MODIS data

We retrieved three snow physical parameters: the snow grain size in the top layer and alower layer as well as the snow impurity concentration. MODIS data were used to test
our algorithms. We found the 1.64 µm band to be sensitive to altitude and humidity, so
we used the 1.24 µm band instead of the 1.64 µm band in the snow parameter retrievals.
Since our retrieved snow parameters are unique and not available as standard products from
MODIS or other similar sensors, a direct comparison is not possible. Availability of groundbased measurement data of snow grain sizes and impurity concentrations for validation
purposes is very limited too. However, we may use the retrieved snow grain size and impurity
concentrations to compute the broadband albedo, which is generally available from satellite
data as well as from ground-based measurements. Because our albedo is calculated from theretrieved snow parameters, the albedo validation can be regarded as an indirect validationof the retrieved snow parameters. Figure 33 shows retrieved snow and aerosol parameters
for day 170, 2005. Corresponding results derived from Aqua MODIS data are very similarand consistent with those derived from Terra MODIS data (see Fig. 38).
In year 2014, we made some improvements in the snow parameter retrieval algorithms.
First, we added the elevation correction to the satellite measured radiances. This correction
will improve the retrieval over high altitude areas, such as Greenland and Antarctica. Second,
we changed the 1.64 µm (MODIS 1.63 µm) channels to 1.05 µm (MODIS 1.24 µm) channel
for snow retrieval, because the 1.64 µm channel is too sensitive to changes in elevation andhumidity. Third, this year we have greatly improved the cloud mask and surface classification
algorithm. This improvement is also a big benefit to snow parameter retrievals, since cloudSGLI Algorithm Description

55

3.2 Application and validation of the C2 algorithm

contaminated pixels and sea-ice pixels will cause unreasonable snow retrievals. The quality ofcloud mask and surface classification is playing a significant role in the snow/ice parameters
retrieval.
In order to quantify these improvements, we re-examined a MODIS image over Greenlandobtained in 2005. Using snow direct albedo, we assessed the quality of retrieved snow
parameters. Figure 33 is an example for MODIS Terra data on day 170 of 2005. In Fig. 34
we compared the retrieved indirect albedo and direct albedo from 2 different versions of ouralgorithm. The top row shows results obtained from the version available in March 2014, andbottom row corresponding results obtained from the improved version of June 2014. The left
column shows the direct albedo, the middle column shows the direct albedo, and the rightcolumn shows the absolute relative difference between the 2 albedo values. In the March2014 version, the cloud mask had some improvements, but the surface classification was not
updated yet. We can see that at the edge of the Greenland ice sheet, there are some bare
ice pixels with low values of the direct albedo, not reflected in the indirect albedo, because
the algorithm falsely retrieved them as “snow pixels”, which led to incorrect albedo values.
By June of 2014, we had upgraded the cloud mask and surface classification algorithm,
and the elevation correction had also been included in snow retrieval algorithm. In thiscurrent version, the indirect albedo and the direct albedo results are very close. Hence, fromthis direct albedo comparison, we infer that there is a significant improvement in our snow
parameter retrievals. Figure 34 is just one example. Retrievals obtained for the whole year
show similar results.

SGLI Algorithm Description

56

3.2 Application and validation of the C2 algorithm

Figure 34: Snow parameter retrievals: Quality check for Terra MODIS images over Greenland, Day 170, 2005 by using different versions of cloud mask, surface classification, andsnow retrieval algorithms. Top: for the version of March 2014; Bottom: for the improved
version of June 2014. Left: indirect albedo; Middle: direct albedo; Right: absolute relative
difference (%).
In year 2015, we updated our snow model by using the non-spherical Voronoi particle
model instead of the spherical particle model. We did a comparison of modeled BRDFs
with NASA aircraft CAR measurements over Barrow, Alaska in 2008 (see Fig. 14). This
comparison showed that the BRDF produced by the use of a Voronoi particle model ismuch closer to the measured snow BRDF than the one obtained from the use of a sphericalparticle model. There is a significant improvement in the phase function shape of Voronoi
particles, which lacks the rainbow patterns produced by a spherical particle model. In order
to correctly use the Voronoi particle phase function, we used a delta-fit technique to get asmooth phase function (see section 2.2.1 and 2.4). In Figs. 35 and 36, we show 2 cases ofsnow retrievals, one based on a spherical snow model and the other on a non-spherical snow
model. These figures clearly show the big improvement obtained from use a Voronoi particle
model. The number of retrievable pixels increased by more than 20%. The image looks moresmooth and the retrieved aerosol values are more reasonable.

SGLI Algorithm Description

57

3.2 Application and validation of the C2 algorithm

Grain size Rs1 (um)(sphere)
Pixel NO.:797258

Grain size Rs2 (um) (sphere)

1000

1000

100

100

Soot Impurity (ppmw)(sphere)
1

Aerosol tau (sphere)
0.05
0.04

0.1

0.03
0.01
0.02
0.001
10

10

0.01

0.00010Grain size Rs1 (um)(nonsphere)
Pixel NO.:1007683 ï 21% more Grain size Rs2 (um) (nonsphere) Soot Impurity (ppmw)(nonsphere)
1
1000

1000

Aerosol tau (nonsphere)
0.05
0.04

0.1

0.03
0.01
100

100

0.02
0.001

10

10

5

4

0.0001

4

x 10

15

0.01
0

4

x 10

3

4

x 10

3

x 10

sphere3

nonsphere

10

2

2

5

1

1

2
1
0

2

0

10

0
ï4
10

2

10

ï2

10

0

10

0

0

0.05

0.1

Figure 35: Snow parameter retrieval over Greenland from a spherical particle model (toppanels), and from a Voronoi particle model (middle panels). The bottom row shows comparison of the retrieved number of pixels.

SGLI Algorithm Description

583.2 Application and validation of the C2 algorithm

Grain size Rs1 (um)(sphere)
Pixel NO.:624008

Grain size Rs2 (um) (sphere)

1000

Soot Impurity (ppmw)(sphere)
1

1000

Aerosol tau (sphere)
0.05
0.04

0.1

0.03
0.01
100

100

0.02
0.001

10

10

Grain size Rs1 (um)(nonsphere)
Pixel NO.:826321 ï 24% more

Grain size Rs2 (um) (nonsphere)

1000

0.01

0.00010

Soot Impurity (ppmw)(nonsphere)
1

1000

Aerosol tau (nonsphere)
0.05
0.04

0.1

0.03
0.01
100

100

0.02
0.001

10

10

5

3

0.01

0.0001

0

4

x 10

4
sphere

2

4

x 10

15000

3

4
3

10000

nonsphere
2

1

0

2
5000

1

2

0

10

x 10

0
ï4
10

2

10

1

ï2

10

0

10

0

0

0.05

0.1

Figure 36: Saimilar retrieval as in Fig. 35 for a MODIS image over Greenland on day 191 of
2008.

SGLI Algorithm Description

59

1
1
−25
0.95
0.95
0.95

0.95

Albedo − SW
Albedo−−SW
SW
Albedo0.95−300.90.9
0.9
−35
0.85
0.85
0.85
−40
0.80.8
100
0.8

120

0.75
0.75
0.75

140

160

180

200

220

240

260

280

Julian day of 2005

−30
0.9
0.9
0.9
−35
0.85
0.85
0.85
−40
0.8
0.8
0.8 100

Albedo − SW
AlbedoAlbedo−−SW
SW

1
1
−25
0.95

−25

−25

−30

−30

−35

120

140

160

180

200

220

240

260

−35

−40

280

100

120

140

160

Julian day of 2005

0.75
0.75
0.75

180

200

220

240

260

−40

280

100

120

140

160

Julian day of 2005

180

Julian day

3.2 Application and validation of the C2 algorithm

−25
−25
−25
0.8
0.8
0.8

−30

−30
−30
−30
0.75
0.75
0.75

−35
0.7−35−40−400.65−400.6
0.6
0.6 100
100
100

0.6

73
176
176

0.6
0.60.6
100
100
100

250
250
250

SummitSummitNASA−SE
SaddleSaddleAlbedo − SW

SurfaceTemp(C)
Albedo− SW
SurfaceTemp(C)
Albedo− −SW
SurfaceTemp(C)
Albedo
SW

250
250
250

150
200
150
200
150
Julianday of200
2005−35−35−35

Albedo− SW
SurfaceTemp(C)
SurfaceTemp(C)
Albedo
− SW

4141
0.8

55

184
0.8 184−15

Temp − AWS
Temp − SIT

55
189
189
0.8

−5−5

0.70.7−10−10150150200200Julianday
2005
Julian
day
of of
2005

250
250

−25
−25
−30
−30
−35
−35

150
150

200
200

Julianday
dayofof2005
2005
Julian

246
246

−15
−15

0.6
0.6

250
250

100
−20100

150
150

−20

−25
−25

−25
−25

−30
−30

−30
−30

200
200

Julianday
dayofof2005
2005
Julian

250
250

−35
−35

−40
−40
100 120 140 160 180 200 220 240 260 280
100 120 140 160 180 200 220 240 260 280

Julian day of 2005
Julian day of 2005

−40
−40 100 120 140
100 120 140

200

160
160

180
180

200
200

220
220

240
240

260
260

280
280

−40
−40 100
100

Surface Temp(C)

250

120
120

140
140

160
160

0

00

−5

0.7

−5−5

−10

−10
−10

0.75

150

−20

−20
−20

Julian day of 2005

−25

−25
−25

−30

−30
−30

0.95

−35

−35
−35

200

−40

−40

5

0.8

0

0.75

JulianJulian
day ofday2005
of 2005
Julian day of 2005

South−D

41

100

150

−20

200

Julian day of 20

100

150

200

−40

0.9

−10
−15
−20

200

220

Julian day of 2005

240

260

280

180

200

NASA−E

10
5

0.65

184

0

0.6
−5

100

150

200

Julian day of 20

−10
−15
−20

−30

180

160

0.7

−25

160

140

0.8

−30

140

120

0.75

−25

−35

100

Julian day of 20

250

Julian day of 2005

120

NASA−SE

0.85

Surface Temp(C)Surface Temp(C)

−15

1

87

0

60

−10

0.6

−35

0.7

SGLI Algorithm Description

0.7

0.95

10

100

284

−30

0.8

−40

Tunu−N10140 140160 160
180 180
200 200
220 220240 260120
240 280
260 280
120 140 160 180 200 220 240 260 280

0.75

−5

SaddleJulianday of20
Julian
day

0.85

NGRIP

100

0.6

200

−25

0.85

0.65

Julian day of 20−40 −40
100 100
120 120140 140160 160
180 180
200

0.65

250

Snow/ice “particle” size and impurity concentration products are new products for EOS.
One could use field measurements to validate the retrieved snow/ice “particle” size andimpurity concentration [Aoki et al., 2007b], but the sparsity of field data, which may beavailable only for a few days of the year at a few locations, represents a significant limitation. However, field albedo measurements are relatively abundant and data sources like theGreenland Climate network (GC-NET) provide important surface measurements including
air temperature, wind speed, wind direction, humidity, pressure, snow depth and albedo
from 18 automatic weather stations (AWS) in Greenland. We will compare the field measured albedo with the retrieved albedo by our snow/ice albedo product. This comparison
will serve as an indirect validation of our snow/ice physical parameter products such as the“particle” size and impurity concentration. Figure 37 shows some preliminary validation
results with GC-NET AWS data for surface short-wave broadband albedo and temperatureduring year 2005. MODIS Terra data have been used by our algorithm (referred to as the
SIT algorithm hereafter) to derive the surface albedo. The results show general agreementof the retrieved surface albedo and temperature with the field measurement.
Another limitation of the field measurement data is that they are obtained at one particular location and do not allow for validation of snow/ice parameter retrievals over large
areas. In view of this space-time sparsity of field measurements and the critical need forvalidations, we propose to use the MODIS/VIIRS broadband albedo products to achieve
this goal. Figure 38 shows a preliminary comparison between our albedo product with the
MODIS MCD43B3 1 km albedo product using MODIS images between June 2 and 10, 2005
over Greenland. In this comparison, we configured our SIT algorithm to retrieve albedo only
5150

−35 −35

−5

115

−15
−15

100

−40 100120
0.9 100

180 200 220 240 260 280
180 200 220 240 260 280

Julian day of 2005
Julian day of 2005

−30 −30
0.6
100

0.9

55
0.8 144 9191

0.6

−20 −20
0.7

1

5

−15

−10 −100.8−15 −150.75

NASA−SE
SummitNASA−SE

1010

10

0.65

200

0.95

0.85

1

−35
−35

Julian day of 2005
Julian day of 2005

0
−5

Julian day of 2005

Snow albedo and temperature – comparisons with field measurement data
Albedo − SW

3.2.3

−10
−10

0.65
0.65

1717

−15
−15

−20

00

−5
−5

0.7
0.7

−10
−10

0.60.6100−20100

0.8

0.75
0.75

0.65
0.65

9292

−15
−15

−20

00

−5
−5

0.70.70.65
0.65

0.60.6100−20 100

0.8

0.75
0.75

Surface Temp(C)

00

0.90

JAR−1
JulianJulian
day
ofday
2005
Juliandayof2005
2005
Julianday
ofof2005

Surface Temp(C)Surface Temp(C)

0.8

180

NASA−
Craeford−1
Summit

10

−5
0.85

304
150

Figure 37: Snow broadband albedo and surface temperature validation against GC-NETAWS data for year 2005 in Greenland. From left to right: Summit, Tunu-N, NASA-E,
Swiss-camp.
0.75
0.75

160

Julian day−25 −25
0.65

0.9

NGRIP
NGRIP

10
10

Albedo − SW

Albedo − SW

Surface
SurfaceTemp(C)Temp(C) Albedo − SW

Surface Temp(C) Albedo − SW
Surface Temp(C) Albedo − SW

55

0.9
0.9

0 00

−40
−40
−40 −40
100
120
140
160
180
200
220
240
260
280
100120120140140160160180180200200220220240240260260280280
100 100
120
140
160
180
200
220
240
260
280

1

0.85
0.85

140

150

5

100

Julianday
Juliandayof2005
2005
PetermanEL
Julian
day
of of
2005
PetermanEL

NASA−E
NASA−E120

0.955 115 91−35−35−35 −35

0.95

10
10

10
1

5

−30
−30
−30
0.6−30

250
250
250

Julianday
dayofof2005
2005
Julian11

0.90.9

SaddleNGRIP
SaddleSwiss−Camp
Humboldt

−20
−20
−20
0.7−20

0.95
0.95

0.85
0.85

100

100

Julian day of 20

−10
−10
−10
0.8−101 1South−DSouth−D

0.6

−25
−25
−25 −25
0.65

246
235
235

−40
−40−40
180
200
220
240
260
280
100120120
120140140
140160160
160180180200200220220240240260260280280
100100Juliandayof 2005
2005
Julianday
JAR−2
Julian
day
of of
2005
JAR−2

250

−15
−15 −15
0.75

0.95
0.95

1010

Julianof 2005
lian
daydayof 2005−25−25−25

Julianday
dayofof2005
2005
Julian

99

200

−5

0.95
0.95

0.90.9250250

150
200
150
150
Julian day 200
of200
2005

Juliandayof2005
2005
Julianday
ofof2005
Julianday

150−5 −5−5
0.85

1 1

0.85
0.85

200200

0.7
0.7

0.65
0.65
0.65
0.6
−30
−30−300.6 100
0.6
100
100

17
4848

−40
−40
−40
100120120
120140140
140160160
160180180
180200200
200220220
220240240260260280280
100100Juliandayofof2005
2005
Julianday
NASA−SE
NASA−SE

87 87

0.7
−20
−20
−20

0.70.7100 120120 140140 160
160 180
180 200
200 220220 240240 260
260 280
280
100

0
0.9

0.8
0.8

0.65
−25
−25
−25

−35
−35
−35

100

5

0.75
−15
−15
−15
0.75
0.75

0.7
−20
−20
−20

250
250
250

0.65 −40

5
5
0.95
6666
112 189

−10
0.8
−10
−10

0.80.8

Julianday
2005
Julian
day
of of
2005

0.7 −35−40
100 120140140160160180 180200 200220 220240 240260 260280 280
100
100 120120 140 160 180 200 220 240 260 280

10101 1010

0 00
0.9
0.9
0.9

0.75
0.75

150
200
150150Julian day
200200of 2005

NASA−E
DYE−2
PetermanEL
NASA−E
NASA−E

−5
0.85
−5 −5
0.85
0.85

0.8 −250.75 −30

Julian day of 2005

5 55 52
0.95
0.95
0.95 18452

0.75
−15
−15
−15

0.65
0.65
0.6
−30
−30
−30
0.60.6 100
100
100

235
67
67

250
250
250

11

−10
0.8
−10
−10

92
7373

200
200
200

JulianJulianday
ofof2005
2005
Julianday
dayof
2005

10
10110

0 0
0
0.9
0.90.90.80.8

Juliandayof2005
2005
Julian
day
of of
2005
Julian
day

150
150
150

0.85 −20

−35

DYE−2

South−DSouth−D

0.85
−5 −5
−5
0.85
0.85

−10
−10 0.8

NGRIP
NGRIP

Tunu−NSouth−D
Tunu−NJAR−2

55
5 0.9541 284
284
0.95
0.95

−15
−150.75
0.75
0.75

−40
−40

250
250
250

1 1

0 0 0.9
0.90.9

−35
−35

200
200
200

Juliandayof2005
Julian
day
ofof
2005
Julian
day
2005

10 10
101

−5−50.85
0.85
0.85

−20
−20 0.7
0.70.7−25−250.65
0.65
0.65
−30
−30 0.6 100
0.60.610010048
330
330

Julianday
Julian
day
of of
2005
Juliandayof 2005
2005

150
150
150

Albedo− SW
SurfaceTemp(C)200
200200Surface Temp(C)Surface Temp(C)Surface Temp(C)JulianJulianday
dayofof2005
2005

Juliandayof2005
2005
Julian
day
of of
2005
Julian
day

55
0.95
144
0.95
0.95 144

Julianof 2005
lian
daydayof 2005

−30

−35
−35
−35
0.7
0.7
0.7

1 1 1

150
200
250
150 Julian 200
200
day
of 2005 250250

−25

−25
0.8−25−30
0.75−30−40−40 100 120 140 160 180 200 220 240 260 280
0.65
0.65
0.65
100100 120120 140140 160160 180
180 200
200 220220 240240 260
260 280280

1501501500.9 −15

−20

−40
−40
−40
0.65100 120 140 160 180 200 220 240 260 280
0.65
0.65
100100120120140140160160180180200200220220240240260260280280

1010

87
115115

−15

−35
−35
−35
0.7
0.7
0.7

Albedo− SW
AlbedoSW
SurfaceTemp(C)
Albedo− −SW
Surface
Surface Temp(C)Temp(C)

NGRIP
JAR−1
JAR−1

−30

Craeford−1−50.95 −10−200.85−20−40

0.60.60.6
100100100

1

SIT
AWS

−10

−15
0.9−15

−25
−25
−25
0.8
0.8
0.8

66

0

Swiss−Camp−100.95−10

−30
−30
0.75
0.75
0.75

Albedo − SW
Albedo − SW

304304115
250250250

5

−5 −5
1 −5

−35
−35
0.7
0.70.7
−40
0.65
0.65
100 120 140 160 180 200 220 240 260 280
0.65
100 120 140 160 180 200 220 240 260 280

Sadd
10

5 5 5 237
52237
0 00

DYE−2
DYE−2
NASA−E

−5
−5 −5
111
−10
−10−10
0.95
0.95
0.95
−15
−15
−15
0.9
0.9
0.9
−20
−20
−20
0.85
0.85
0.85

Humboldt
DYE−2
Humboldt

101010

Albedo − SW

−25
−25
0.8
0.80.8

−30
0.75
0.75
0.75

182
67182
250
250
250

Albedo − SW

Temp − AWS
Temp − AWS
Temp − SIT
Temp − SIT

−20

200
200
200

Julianday
ofof2005
Julianday
dayof
2005
Julian2005

Surface Temp(C)

−15
−15
0.9
0.90.9

150
150
150

5 55 143
2841430 00

Tunu−N
Tunu−NSouth−D

−5 −5
1−511
−10
−10
−10
0.95
0.95
0.95
−15
−15
−15
0.9
0.9
0.9
−20
−20
−20
0.85
0.85
0.85

NASA−U
Tunu−NNASA−U

10

10 10
0.6
0.6
0.6
100
100
100

250
250
250

Surface Temp(C)

1

−20
0.85
0.85
0.85

200200200

244
244
330

200
200
200

Julianday ofof2005
JulianJulian day
day of 2005
2005

Albedo − SW

Albedo − Surface
SW
Temp(C)
Albedo −Surface
SW Temp(C)Albedo − SW

SummitSummitSaddle−10

Julianof 2005
lian
daydayofday
2005
Julianof 2005

150
150
150

0 00

0

−5
1−5 1

−10
0.95
0.95
0.95

150 15010100.60.6
10
0.6
100
100
100
5 55 115
144115

250
250
250

Surface Temp(C)

132
132
176

200200200

Juliandayof 2005
Julianday
Juliandayofof2005
2005

Surface
SurfaceTemp(C)Temp(C)0

Humboldt
Humboldt
JAR−1

150150150

Albedo − SW

5 112
112

0.65
0.65
0.65

Albedo −−Surface
SW
AlbedoSW
Temp(C)
Albedo− SW

5Craeford−1
SummitCraeford−1

Surface Temp(C)Surface Temp(C)Surface Temp(C)

10

0.610
0.60.6100100
100

0.7
0.7
0.7

0.65
0.65
0.65

SurfaceTemp(C)Temp(C)
Surface304

250250250

0.70.7
0.7

Swiss−Camp
Swiss−Camp

Temp(C)
Albedo −Surface
SW
SurfaceTemp(C)Temp(C)
AlbedoSW
Albedo −−SurfaceSW

99 99200200200

liandaydayof 2005
Julianof 2005
Juliandayof 2005

0.65
0.65
0.65

Albedo − SW

0.70.70.7

150 150

150

1
1
0.95
0.95
0.95
0.90.9
0.9
0.85
0.85
0.85
0.80.8
0.8
0.75
0.75
0.75

Albedo − SW
Albedo − SW
Albedo − SW

SITSIT
AWS
AWS

−35
−40

100

120

140

160

180

200

Julian day of 20

3.2 Application and validation of the C2 algorithm

for snow pixels, so that land and ice pixels on the melted and melting part of Greenland are
not included in our results. Figure 38 shows that there is basic agreement and we will refine
it into a pixel by pixel comparison for quantitative and statistical validation.

Figure 38: Albedo retrieval compared to MODIS MCD43 surface albedo between June 2
and 10, 2005 over Greenland. Upper panels: MCD43; Lower panels: SIT algorithm. Fromleft to right: retrieved albedo for VIS (0.3-0.7 µm), NIR (0.7-2.8 µm) and SW (0.3-2.8 µm)
wavelength ranges.

3.2.4

Summary

In JFY 2015, we completed the 3-year project. There were several improvements, upgrades
and tests. The most important results obtained this year may be summarized as follows:
1. We acheived a very significant improvement in the cloud mask. The upgraded algorithm has been validated by use of CALIOP data, and shown to have better performance than the MODIS algorithm.
2. The surface classification algorithm has also been upgraded, especially for the sea
ice/water and snow-covered sea ice/bare sea ice separation. The dynamic thresholdmethod will also improve the accuracy of classification under different solar/viewing
geometries.
3. A non-spherical Voronoi particle snow model has been implemented in the snow retrieval algorithm. A comparison of snow BRDFs obtained from model simulations
and NASA CAR measurements shows that the BRDF produced by a Voronoi particle
SGLI Algorithm Description

613.2 Application and validation of the C2 algorithm

model is closer to the measured snow BRDF than the one produced by a spherical particle model. Further, MODIS-retrieved snow parameters based on a Voronoi particle
model are more accurate than those based on a spherical particle model.
4. Using response functions for the SGLI channels, we generated snow/ice surface temperature retrieval tables for the SGLI sensor, as we did for MODIS. The new tables
have been implemented in the SGLI snow retrieval code.
5. We did a sensitivity study to explore the possibility of improving aerosol retrievals
over snow by using the SGLI polarization channels. The results show that when the
scattering angle is less than 90 degrees, the polarized reflectance could be used to
distinguish effects of aerosols from those of snow impurity. The SGLI polarized sensor
orbit direction design satisfy this condition. So use of polarized SGLI measuremenst
could lead to improved aerosol retrievals over snow. We will continue this work in thenext 3-year project.
6. In the next 3-year RA6 project, we will continue to study the non-spherical particle
snow model and compare results with field measurements. Also we will explore thebenefits obtained by using SGLI polarized channels and multiple angle measurements
in snow retrievals. After SGLI launch, we will work on testing of the snow/ice products
and further validations.

SGLI Algorithm Description

62

REFERENCES

References
Ackerman, S. A., Strabala, K. I., Menzel, W. P., Frey, R. A., Moeller, C. C., & Gumley,
L. E. (1998). Discriminating clear sky from clouds with MODIS. Journal of GeophysicalResearch, 103 , 32141.
Anderson, G. P., Clough, S. A., Kneizys, F. X., Chetwynd, J. H., & Shettle., E. P. (1986).
AFGL atmospheric constituent profiles (0-120km), AFGL-TR-86-0110 (OPI). Optical
Physics Division, Air Force Geophysics Laboratory Hanscom AFB, MA 01736.
Aoki, T., Aoki, T., Fukabori, M., Hachikubo, A., Tachibana, Y., & Nishio, F. (2000). Effects
of snow physical parameters on spectral albedo and bi-directional reflectance of snow
surface. Journal of Geophysical Research, 105 , 10,219 – 10,236.
Aoki, T., Hori, M., Motoyohi, H., Tanikawa, T., Hachikubo, A., Sugiura, K., Yasunari, T. J.,
Storvold, R., Eide, H. A., Stamnes, K., Li, W., Nieke, J., Nakajima, Y., & Takahashi, F.
(2007a). ADEOS-II/GLI snow/ice products: Part II - Validation results using GLI and
MODIS data. Remote Sensing of Environment, 111 , 274 – 290.
Aoki, T., Hori, M., Motoyoshi, H., Tanikawa, T., Hachikubo, A., Sugiura, K., Yasunari, T. J.,
Storvold, R., Eide, H. a., Stamnes, K., Li, W., Nieke, J., Nakajima, Y., & Takahashi, F.
(2007b). ADEOS-II/GLI snow/ice products Part II: Validation results using GLI and
MODIS data. Remote Sensing of Environment, 111 , 274–290.
Barton, I. J. (1985). Transmission model and ground-truth investigation of satellite-derived
sea surface temperatures. J. Clim. Appl. Meteorol., 24 , 508 – 516.
Blanchet, J. P., & List, R. (1983). Estimation of optical properties of arctic haze using a
numerical model. Atmos. Ocean., 21 , 444–465.
Chen, N., Li, W., Tanikawa, T., Hori, M., Aoki, T., & Stamnes, K. (2014). Cloud Mask
Over Snow/ice Covered Areas for the GCOM-C1/SGLI Cryosphere Mission: Validations
over Greenland. Journal of Geophysical Research: Atmospheres, 119 .
Dahlback, A., & Stamnes, K. (1991). A new spherical model for computing the radiation
field available for photolysis and heating at twilight. Planetary and Space Science, 39 ,671–683.
Dozier, J., & Warren, S. G. (1982). Effect of viewing angle on the infrared brightness
temperature of snow. Water Resour. Res., 18 , 1424–1434.
Ettema, J., van den Broeke, M. R., van Meijgaard, E., van de Berg, W. J., Box, J. E., &Steffen, K. (2010). Climate of the Greenland ice sheet using a high-resolution climate
model – Part 1: Evaluation. The Cryosphere, 4 , 511–527.
Gjerstad, K. I., Stamnes, J. J., Hamre, B., Lotsberg, J. K., Yan, B., & Stamnes, K. (2003).
Monte carlo and discrete-ordinate simulations of irradiances in the coupled atmosphereocean system. Applied Optics, 42 , 2609 – 2622.
SGLI Algorithm Description

63

REFERENCES

Grenfell, T. C., & Warren, S. G. (1999a). Relationship of a nonspherical ice particle by acollection of independent sphere for scattering and absorption of radiation. Journal ofGeophysical Research, 104 , 31,697 – 31,709.
Grenfell, T. C., & Warren, S. G. (1999b). Representation of a nonspherical ice particle bya collection of independent sphere for scattering and absorption of radiation. Journal ofGeophysical Research, 104 , 31697–31709.
Hall, D. K., Box, J. E., Casey, K. A., Hook, S. J., Shuman, C. A., & Steffen, K. (2008). Comparison of satellite-derived and in-situ observations of ice and snow surface temperatures
over greenland. Remote Sensing of Environment, 112 , 3739–3749.
Hall, D. K., Riggs, G. A., & Barton, J. S. (2001). Algorithm Theoretical Basis Document( ATBD ) for the MODIS Snow and Sea Ice-Mapping Algorithms. Technical Report
September.
Hall, D. K., Riggs, G. A., & Salomonson, V. V. (1995). Development of Methods for Mapping
Global Snow Cover Using Moderate Resolution Imaging Spectroradiometer Data. Remote
Sensing of Environment, 34 , 127–140.
Hall, D. K., Riggs, G. A., & Salomonson, V. V. (1996). Algorithm Theoretical Basis Document ( ATBD ) for the MODIS Snow and Sea Ice-Mapping Algorithms. Technical Report
September.
Hamre, B., Winther, J.-G., Gerland, S., Stamnes, J. J., & Stamnes, K. (2004). Modeled
and measured optical transmittance of snow-covered first-year sea ice in kongsfjorden,svalbard. Journal of Geophysical Research: Oceans, 109 .
Hanssen, A., & Kuipers, W. (1965). On the Relationship Between the Frequency of Rainand Various Meteorological Parameters: (with Reference to the Problem Ob Objective
Forecasting). Mededelingen en verhandelingen. Staatsdrukkerij- en Uitgeverijbedrijf.
Hess, M., Koepke, P., & Schult, I. (1998). Optical properties of aerosols and clouds: The
software package opac. Bulletin of the American Meteorological Society, 79 , 831 – 844.
Holz, R. E., Ackerman, S. A., Nagle, F. W., Frey, R., Dutcher, S., Kuehn, R. E., Vaughan,
M. A., & Baum, B. (2008). Global Moderate Resolution Imaging Spectroradiometer(MODIS) cloud detection and height evaluation using CALIOP. Journal of GeophysicalResearch, 113 , D00A19.
Hori, M., Aoki, T., Tanikawa, T., Motoyoshi, H., Hachikubo, A., Suqiura, K., T J Yasunari,
H. E., Storvold, R., Nakajima, Y., & Takahashi, F. (2006). In-situ measured spectral
directional emissivity of snow and ice in the 8?14 µm atmospheric window. Remote Sensing
of Environment, 100 , 486?502.
Hu, Y. X., & Stamnes, K. (1993). An Accurate Parameterization of the Radiative Properties
of Water Clouds Suitable for Use in Climate Models. Journal of Climate, 6 , 728–742.

SGLI Algorithm Description

64

REFERENCES

Hutchison, K. D., Iisager, B. D., & Mahoney, R. L. (2013). Enhanced snow and ice identification with the VIIRS cloud mask algorithm. Remote Sensing Letters, 4 , 929–936.
Jiang, S., Stamnes, K., Li, W., & Hamre, B. (2005). Enhanced solar irradiance across theatmosphere–sea ice interface: a quantitative numerical study. Appl. Opt., 44 , 2613–2625.
Jin, Z., Stamnes, K., Weeks, W. F., & Tsay, S.-C. (1994). The effect of sea ice on thesolar energy budget in the atmosphere-sea ice-ocean system: A model study. Journal ofGeophysical Research: Oceans, 99 , 25281–25294.
Key, J. R., Collins, J. B., Fowler, C., & Stone, R. S. (1997). High-latitude surface temperature
estimates from thermal satellite data. Remote Sens. Environ., 67 , 302 – 309.
Key, J. R., & Haefliger, M. (1992). Tba. J. Geophys. Res., 67 , 5885 – 5893.
Key, J. R., Yang, P., Baum, B. A., & Nasiri, S. L. (2002). Parameterization of shortwave
ice cloud optical properties for various particle habits. Journal of Geophysical Research,
107 , 1–10.
Kuchiki, K., Aoki, T., Niwano, M., Motoyoshi, H., & Iwabuchi, H. (2011). Effect of sastrugi on snow bidirectional reflectance and its application to MODIS data. Journal ofGeophysical Research, 116 , D18110.
Li, W., Stamnes, K., Eide, H., & Spurr, R. (2007). Bidirectional reflectance distributionfunction of snow: corrections for the lambertian assumption in remote sensing applications.
Optical Engineering, 46 , 066201–066201–9.
Li, W., Stamnes, K., Spurr, R., & Stamnes, J. (2008). Simultaneous retrieval of aerosol
and ocean properties by optimal estimation: Seawifs case studies for the santa barbara
channel. International Journal of Remote Sensing, 29 , 5689 – 5698.
Liu, Y., Ackerman, S. A., Maddux, B. C., Key, J. R., & Frey, R. A. (2010). Errors inCloud Detection over the Arctic Using a Satellite Imager and Implications for Observing
Feedback Mechanisms. Journal of Climate, 23 , 1894–1907.
Llewellyn-Jones, D. T., Minnett, P. J., Saunders, R. W., & Zavody, A. M. (1984). Satellite
multichannel infrared measurements of sea surface temperature of the northeast atlantic
ocean using avhrr/2. Q. J. R. Meteorol. Soc., 110 , 613 – 631.
McClain, E. P., Pichel, W. G., & Walton, C. C. (1985). Comparative performance of avhrrbased multichannel sea surface temperature. J. Geophys. Res., 18 , 11587–11601.
Miller, S. D., & Lee, T. F. (2005). Satellite-Based Imagery Techniques for Daytime Cloud /Snow Delineation from MODIS. Journal of Applied Meteorology, 44 , 987–997.
Minnett, P. J. (1990). The regional optimization of infrared measurements of sea surface
temperature from space. J. Geophys. Res., 95 , 13497 – 13510.

SGLI Algorithm Description

65

REFERENCES

Price, J. C. (1983). Estimating surface temperature from satellite thermal infrared data – asimple formulation for the atmospheric effect. Remote Sens. of Environ., 13 , 353–361.
Rodgers, C. D. (2000). Inverse methods for atmospheric sounding - Theory and Practice -.
Series on Atmospheric, Oceanic and Planetary Physics. World Scientific.
Salomonson, V., & Appel, I. (2004). Estimating fractional snow cover from MODIS using
the normalized difference snow index. Remote Sensing of Environment, 89 , 351–360.
Salomonson, V., & Appel, I. (2006). Development of the Aqua MODIS NDSI fractional snow
cover algorithm and validation results. IEEE Transactions on Geoscience and Remote
Sensing, 44 , 1747–1756.
Spurr, R. (2002). Simultaneous derivation of intensities and weighting functions in a general
pseudo-spherical discrete ordinate radiative transfer treatment. Journal of Quantitative
Spectroscopy and Radiative Transfer , 75 , 129 – 175.
Stamnes, K., Hamre, B., Stamnes, J., Ryzhikov, G., Biryulina, M., Mahoney, R., Hauss, B.,& Sei, a. (2011). Modeling of radiation transport in coupled atmosphere-snow-ice-oceansystems. Journal of Quantitative Spectroscopy and Radiative Transfer , 112 , 714–726.
Stamnes, K., Li, W., Eide, H., Aoki, T., Hori, M., & Storvold, R. (2007a). ADEOS-II/GLI
snow/ice products Part I: Scientific basis. Remote Sensing of Environment, 111 , 258–273.
Stamnes, K., Tsay, S. C., Wiscombe, W., & Jayaweera, K. (1988a). Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting
layered media. Applied optics, 27 , 2502–9.
Stamnes, K., Tsay, S. C., Wiscombe, W. J., & Jayaweera, K. (1988b). Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering andemitting layered media. Applied Optics, 27 , 2502 – 2509.
Stamnes, K., Wei, L., Eide, H., Aoki, T., Hori, M., & Storvold, R. (2007b). ADEOS-II/GLI
snow/ice products: Part I - Scientific basis. Remote Sensing of Environment, 111 , 258 –273.
Trepte, Q., Minnis, P., & Arduini, R. F. (2002). Daytime and nighttime polar cloud andsnow identification using MODIS data. In Proc. SPIE 3rd International Asia-Pacific
Environmental Remote Sensing Symposium: Remote Sensing of the Atmosphere, Ocean,
Environment, and Space, Hangzhou, China, October 23-27 (pp. 449–459).
Wan, Z. (2008). New refinements and validation of the modis land-surface temperature/emissivity products. Remote Sensing of Environment, 112 , 59?74.
Wan, Z., & Dozier, J. (1996). A generalized split-window algorithm for retrieving landsurface temperature from space. IEEE Trans. Geosci. Remote Sens., 34 , 892–905.
Warren, S. G. (1982). Optical properties of snow. Reviews of Geophysics and Space Physics,
20 , 67 – 89.
SGLI Algorithm Description

66

REFERENCES

Warren, S. G., & Brandt, R. E. (2008). Optical constants of ice from the ultraviolet tothe microwave: A revised compilation. Journal of Geophysical Research, 113 , D14220.
doi:10.1029/2007JD009744.
Warren, S. G., Brandt, R. E., & Hinton, P. O. (1998). Effect of surface roughness on
bidirectional reflectance of Antarctic snow. Journal of Geophysical Research, 103 , 25789–
25807.

SGLI Algorithm Description

67