NDX-120015A

Descriptions of GCOM-W1 AMSR2Level 1R and Level 2 Algorithms

Japan Aerospace Exploration Agency
Earth Observation Research Center

Jul. 9, 2013

NDX-120015A

Chapter 1. Level 1R
Chapter 2. Level 2 (Integrated Water Vapor and Cloud Liquid Water)
Chapter 3. Level 2 (Precipitation)
Chapter 4. Level 2 (Sea Surface Temperature)
Chapter 5. Level 2 (Sea Surface Wind Speed)
Chapter 6. Level 2 (Sea Ice Concentration)
Chapter 7. Level 2 (Snow Depth)
Chapter 8. Level 2 (Soil Moisture)

NDX-120015A

Chapter 1.
Description of GCOM-W1 AMSR2 Level 1R productTakashi Maeda
Japan Aerospace Exploration Agency
Earth Observation Research Center

Yuji Taniguchi
Mitsubishi Space Software Co., Ltd.

1-1

NDX-120015A

1. Introduction
In AMSR-E, the footprint size of each frequency can differ within a range of several tens ofkilometers. In addition, the center of a footprint at each frequency that corresponds to the same pixel
and scan numbers can differ within a range of several kilometers owing to the difference in thelocations of the horns of the main rotating reflector. In a precise sense, therefore, even the brightness
temperatures (T B ) of each frequency, which are stored in the same pixel and scan numbers of theL1B product, can be measured in a different area. It was considered that these different measurement
conditions decreased the accuracy of a higher level algorithm when retrieving physical values using
the L1B product as an input.
AMSR2 uses a similar measurement mechanism to AMSR-E, and it experiences a decline in thefootprint size and center at different frequencies. To compensate for this decline, we convolve some
footprints (antenna patterns) for one frequency (f 1 ) using the Backus-Gilbert method and calculate
the weighting factors to obtain a corrected footprint for f 1 , which agreed quite well with the footprint
size and center of another frequency (f 2 ). Here, the T B values of f 1 and f 2 that correspond to the same
pixel and scan numbers are denoted by T B1 and T B2 , respectively. We combine the T B1 values^

observed in the footprints of f 1 using the weighting factors to obtain the corrected T B1 (T B1 ) that is^

expected to be observed in the footprint of f 2 . If T B1 is used in place of T B1 in the phase of the higher
level algorithm, where T B1 was used with T B2 as an input, it is expected that its retrieval accuracy^

will be improved because the footprint centers and sizes of T B1 and T B2 are standardized.
The AMSR2 observation frequencies are 6.9, 7.3, 10.65, 18.7, 23.8, 36.5, and 89.0 GHz. The
datasets of the corrected T B are provided as the L1R products for the vertically and horizontally
polarized components of f 1 and f 2 , which are indicated by O and O* in Table 1.
Table 1: Combination of two frequencies where the dataset for the corrected T B is provided as theL1R product. In this table, the T B of f 1 is corrected as observed in the footprint of f 2 .
f 2 [GHz]f 1 [GHz]

6.9 and 7.3

1

1

10.65

23.8

36.5

6.9

O*

X

X

X

7.3

O*

X

X

X

10.65

O

O*

X

X

18.7

O

O

O

X

23.8

O

O

O*

X

36.5

O

O

O

O*

89

O

O

O

O

The footprint location and size at 7.3 GHz was very similar to that at 6.9 GHz.
1-2

NDX-120015A^

In the L1R product, the footprint center of T B1 is always modified to that of T B at 89.0 GHz (A)^that corresponds to the same pixel and scan numbers. Therefore, the spatial resolution of T B1 is nodifferent from that of T B1 in the combinations indicated by O* in Table 1. In the combinations
indicated by O, the footprint size of f 1 is greater than that of f 2 . In particular, the footprint size at
23.8 GHz is actually slightly greater than that at 18.7 GHz.
At frequencies higher than 7.3 GHz (18.7 GHz), the dataset for T B corrected for the footprint sizeand center at 7.3 GHz (18.7 GHz) is not provided in the L1R product. Therefore, when a T B of 7.3
GHz (18.7 GHz) was used with that of another higher frequency, it is recommended that a T B is usedwith a correction for the footprint size and center at 6.9 GHz (23.8 GHz).

2. Theoretical Description
Smaller antenna patterns with a higher frequency G i (ρ) were convolved to a larger antenna
pattern with a lower frequency F(ρ) by using the weighting factors a i (i = 1, …, N). The original T B^

observed in G i (ρ) was T Bi , whereas the corrected T B that was expected to be observed in F(ρ) was T
B . First, G i (ρ) and F(ρ) are normalized as∫∫ G ( ρ )dA = 1(1)∫∫ F ( ρ )dA = 1

(2)

A

i

A

According to the Backus-Gilbert method, the weighting factor a is calculated as follows.

 1 − u T V −1 v  
a = V −1  v +  T −1 u 
 u V u  (3)

V = G +ε I

(4)u i = ∫∫ Gi( ρ )dA ≡ 1(5)

vi = ∫∫ Gi ( ρ ) F ( ρ )dA(6)

Gij = ∫∫ Gi ( ρ )G j ( ρ )dA

(7)

where

A

A

A

I is the identity matrix and ε is a parameter that determines the features of a. a satisfies the equation
below.
NN∫∫ ∑ a G ( ρ )dA = ∑ a = 1
A

i =1

i

i

i =1

^

Finally, T B is calculated as follows.

1-3

i

(8)

NDX-120015A

N

TˆB = ∑ ai TBi

(9)

i =1

The convolution error is evaluated using the value Q 0 , given below.
2N
Q0 = ∫∫ ∑ ai Gi ( ρ ) − F ( ρ ) dA
A i =1(10)

A prior evaluation using the AMSR-E data showed that when ε is less than 10-5, some factors of a^become less than 0, and Q 0 becomes smaller. This means that some T Bi are subtracted to obtain T B ,^so the calculation to obtain T B loses its physical meaning. On the other hand, when ε > 10-3, Q 0
increases. Therefore, ε is always determined as 10-4 in the L1R product.

3. Algorithm Implementation
Provided that the altitude and attitude of the GCOM-W1 satellite do not change, the distributionof footprints at each frequency is fixed during each scan. Therefore, the altitude and attitude of
GCOM-W1 are assumed to be as shown in Fig. 1. The weighting factors a for each pixel number
during one scan are calculated in advance. a is calculated for a specific pixel number during the scanand the same pixel number is used in other scans. a will be recalculated when the assumption shownin Fig. 1 changes substantially.
GCOM-W1/AMSR2
Off nadir angle 47.5°1124 km

Altitude of GCOM-W1
705 kmEarth’s surface

828.7 kmEarth’s radius
6373 km

6373 km7.48°

Fig. 1 Altitude and attitude of the GCOM-W1 satellite is assumed to calculate the weighting factorsa.

1-4

NDX-120015A

The antenna patterns G i (ρ) convolved to obtain F(ρ) for each pixel P(x, y) are derived from thepixels P(x + Δp, y + Δl) (-30 ≤ Δp ≤ 30, -30 ≤ Δl ≤ 30) (Fig. 2). Thus, the number of G i (ρ) convolved
to obtain F(ρ) (N) is 3721 (= 61 × 61). In any combination of G i (ρ) arnd F(ρ), these G i (ρ) aroundF(ρ) almost cover the area where F(ρ) becomes larger than -20 dB of the maximum value of F(ρ)
wherever F(ρ) is located in the swath, so it is expected that the a is calculated appropriately.

61 pixels
Pixel corresponding to F(ρ)

61 scans

about 2000 scans

243 pixels (486 pixels)1

Pixels corresponding to Gi(ρ)
convolved to obtain F(ρ)

Fig. 2 Distribution of G i (ρ) convolved to calculate F(ρ)^

The corrected T B that is expected to be observed in F(ρ) (T B ) is calculated using Eq. (9) with thea and T B values observed in G i (ρ) (T Bi ). If some values of T Bi are recognized as abnormal, or if theycannot be defined because P(x, y) is located at the edge of the coordinate space composed of the^

pixel and scan numbers, the weighting factors a i corresponding to such T Bi are invalidated to 0, so T^

B can be calculated. However, when the number of invalid T Bi exceeds N/2, T B itself is regarded as

an abnormal value.
As described in Section 1, the latitude λ and longitude φ of the footprint's center at eachfrequency corresponding to P(x, y) are always modified to those at 89.0 GHz (A). However, during^

the actual final phase used to create the L1R product, the geometric location of each T B is correctedbased on the altitude H that corresponds to the λ and φ. To determine H, a dedicated altitude databaseis constructed from GTOPO30 with 0.02° and 0.02° intervals in latitude and longitude, respectively,in advance, whereas H is calculated by averaging the altitude values stored in the database for thearea of 50 km in a radius centered on (λ, φ).

1-5

NDX-120015A

Using this H value, the corrected geometric location approaches to AMSR2 according to

L = H tan θ(11)where θ is the earth incident angle of AMSR2 at the original latitude and longitude (Fig. 3).
GCOM-W1 / AMSR2Corrected latitude and longitude

L

H

θ

Ground level

Original latitude λ and longitude φEarth’s ellipsoid heightFig. 3 Altitude correction during the final phase

References[1] Backus, G. and F. Gilbert, Uniqueness in the inversion of inaccurate gross earth data, Phil. Trans.
Roy. Soc. London, A266, pp. 123-192, 1970.
[2] Stogryn, A., Estimates of brightness temperatures from scanning radiometer data, IEEE Trans. on
Antennas Propagation, AP-26, pp. 720-726, 1978.
[3] Poe, G. A., Optimum interpolation of imaging microwave radiometer data, IEEE Trans. on
Geoscience and Remote Sensing, GE-28, pp. 800-810, 1990.

1-6

NDX-120015A

Chapter 2.
Description of GCOM-W1 AMSR2Integrated Water Vapor and Cloud Liquid Water Retrieval
Algorithm

Masahiro Kazumori
Numerical Prediction Division
Japan Meteorological Agency
2-1

NDX-120015A

Table of Contents

1.

Introduction ................................................................................................................... 3

2.

Algorithm Overview ....................................................................................................... 3

3.

Theoretical Description .................................................................................................. 4

4.

Algorithm Implementation ............................................................................................. 7

5.

4.1.

Implementation ......................................................................................................... 7

4.2.

Input/output parameters ............................................................................................. 7

4.3.

Ancillary data ........................................................................................................... 7

4.4.

Processing flow ........................................................................................................ 8

4.5.

Example output......................................................................................................... 9

4.6.

Limitations ............................................................................................................... 9

Validation Concept ......................................................................................................... 9

References ........................................................................................................................... 102-2

NDX-120015A

1. Introduction
This document describes an algorithm for retrieval of integrated water vapor and cloud liquid waterfrom the Advanced Microwave Scanning Radiometer-2 onboard the GCOM-W1 satellite, which isthe first satellite of Japan’s Global Change Observation Mission (GCOM-W1/AMSR2). Verticallyintegrated atmospheric water vapor and cloud liquid water are key geophysical parameters tounderstanding the global change of the water cycle, and climate change on Earth. Atmospheric water
vapor is the source of atmospheric convection and cloud formation and it plays a crucial role inprecipitation processes. Cloud liquid water is linked to the radiation budget of the Earth, and its
distribution is one of the undetermined parameters in weather and climate prediction modeling for
global warming. Therefore, accurate and long term measurements of these fundamental geophysicalparameters are essential, not only for the environmental monitoring, but also for further development
and verification of numerical weather and climate prediction models. The measurement of theseparameters is one of the main targets of the GCOM mission.
In order to produce accurate and long-term products of these geophysical parameters, an algorithm
is developed based on a physical-statistical model. The algorithm is designed to be applicable forvarious historical microwave radiometers. Moreover, with a view to operational application, thealgorithm is computationally fast for generation of retrieval products in real time. The originalalgorithm (Takeuchi, 2002) that was developed at the Japan Aerospace Exploration Agency (JAXA)
produces a standard AMSR-E integrated water vapor product. The algorithm has been improved anda new capability has been introduced to retrieve cloud liquid water from GCOM-W1/AMSR2.

2. Algorithm Overview
AMSR2 observes Earth radiation from space. Microwave measurements can provide information
about moisture over oceans. Measurements at microwave frequencies, i.e., 19 GHz, 23 GHz, and 37
GHz channels, allow us to observe the emission of atmospheric water vapor and cloud liquid water.
The radiance at the top of the atmosphere contains integrated information on the amount of
atmospheric water vapor and cloud liquid water. Accurate estimation of the emission from the ocean
surface (i.e., estimating ocean emissivity) is important to achieve an accurate retrieval. The ocean
emissivity can be described in terms of microwave frequency, sea surface temperature (SST), sea
surface wind (SSW), and incidence angle. In order to save on the huge calculation cost associatedwith scattering by hydrometeors in the atmosphere, the algorithm focuses on retrieval under rain freeconditions over the ocean. As each AMSR2 channel has a different sensitivity to the water vapor andcloud liquid water, AMSR2’s multi-channel ability and dual polarization measurements are utilized
2-3

NDX-120015A

in the algorithm. The algorithm is capable of retrieving simultaneously the integrated water vapor(Total Precipitable Water, TPW) and Cloud Liquid Water (CLW) over the ocean from the AMSR2
brightness temperature. The calibration of TPW products is based on integrated water vapor fromradiosonde observation (RAOB).

3. Theoretical Description
The algorithm is capable of retrieving TPW and CLW over the ocean from the AMSR2 brightness
temperature under rain free conditions. The retrieval algorithm uses a physical-statistical method andis based on a theoretical microwave radiative transfer model (Stephens, 1994 and Jannsen, 1993).
The model is somewhat simplified in order to be computationally fast for operational applications. A
single-layer atmosphere above the sea surface is considered in the model. A vertically homogeneousand isothermal atmosphere is assumed. Emission and absorption by atmospheric water vapor,oxygen, and ozone, and emission and absorption at the sea surface are simulated in the model. The
effect of varying the salinity of sea water for ocean emissivity (in other words, the dielectric
constant) is not considered because the microwave frequency range of AMSR2 is not very sensitive
to the salinity. The effect of a rough surface ruffled by wind and the variation of ocean emissivity
due to SST are considered in the algorithm.
In general, the radiative transfer equation in the microwave range is written in terms of brightness
temperature as follows:

Tb = Tu + [ε sTs + (1 − ε s ) Td ] Tr

[1]

where ε s is the ocean emissivity, Ts is the SST, and Tr is the atmospheric transmittance. On the righthand side of the equation, the first term denotes upward emission from the atmosphere, the second
term denotes upward emission from the ocean, and the third term denotes the ocean reflection of thedownward emission from the atmosphere.
The mean emission temperature of the atmosphere, Ta , is introduced, and it is assumed that themean temperatures for water vapor and for cloud liquid water in the atmosphere are equal.
Consequently, the mean emission temperature is a parameter in the modeling of the atmosphere thatis not separated into vapor and liquid phases. As a vertically homogeneous and isothermal
atmosphere is assumed, the upward and the downward emission can be set as equal to the meanemission temperature. Under these assumptions, the equation can be written as

Tb = (1 − Tr ) Ta + Tr ε sTs + Tr (1 − ε s ) (1 − Tr ) Ta

[2]

Three different kinds of look-up tables (LUT) are used in the algorithm. The first one is for the
2-4

NDX-120015A

computation of ocean emissivity. The inputs of the LUT are SST, SSW, and the AMSR2 channel
microwave frequencies. The second LUT is for the computation of atmospheric mean emission
temperatures. The inputs of the LUT are the square of the atmospheric transmittance estimated in thealgorithm and the temperature at 850 hPa from Japan Meteorological Agency (JMA) global
numerical weather prediction (JMA 2007) results. The last LUT is for conversion from index
parameters to TPW and CLW.
The ocean emissivity can be obtained from the first LUT by using ancillary SST and SSW data.
The LUT is produced by using a match-up data set of atmospheric profile data from RAOB, AMSR2
brightness temperature observations, and ancillary SST and SSW data. The LUT provides the ocean
emissivity for each AMSR2 frequency as a function of SST and SSW. For this algorithm, the SST
and SSW used as ancillary input data are obtained from JMA global numerical weather prediction
results. For near real time TPW and CLW generations, SST and SSW are obtained from JMA’s
short-range global forecasts. Temperatures at 850 hPa are also obtained from JMA data and are used
as an input in the estimation of atmospheric transmittance at a later stage in the algorithm. Since
TPW and CLW are retrieved for each Field of View (FOV) of AMSR2, these ancillary data wereinterpolated spatially and temporally for the AMSR2 FOV location.
It is relatively easy to describe the characteristics of an ocean surface compared with those of a land
surface. Ocean surface emissivity is characterized by observational microwave frequency, SST, SSW,
incidence angle, and ocean salinity. The dependency of sea surface reflectivity and emissivity on
SSW and SST is not described explicitly in the model. Indeed, the ocean surface is assumed to be a
specular surface and the Fresnel law is used to find the reflectivity and the emissivity of the ocean.
The deviation from the Fresnel law of the rough ocean surface is estimated empirically by using thesecond LUT, produced by matching up data from radiosonde observation, SST and SSW, and
AMSR2 data. The atmospheric profiles from RAOB, SST, SSW and a general microwave radiativetransfer calculation were used to produce the LUT for the empirical emissivity estimation. Thisempirical emissivity estimation method (Kazumori et al., 2008) is based on real observed microwave
brightness temperature. The resulting LUT can capture SST and SSW dependency accurately. Itshould be noted that the accuracy is highly dependent on the calibration accuracy of the input
AMSR2 brightness temperature.
The previous equation can be written by introducing a parameter α as indicated below:

{

}

[3]

(Ts − Ta ) Trε s 2
1 − (1 − ε s ) Tr 

[4]

Tb = α 1 − (1 − ε s ) Trα ≡ Ta +

The square of the transmittance is written as below:

2-52NDX-120015A

Tr =2

1 − Tb / α1− ε s[5]

Figure 1 shows a schematic diagram of the TPW and CLW retrieval algorithm. The transmittance isestimated in the algorithm. In step 1, the ocean emissivity is estimated by using the first LUT,
ancillary SST and SSW data, and AMSR2 channel frequencies at each AMSR2 observation location.
In step 2, an initial value is set for the transmittance. In step 3, the mean emission temperature is
obtained from the second LUT. The transmittance and temperature at 850 hPa at the AMSR2
observation location are used as input data for the LUT. The pre-defined LUT was created inadvance by using atmospheric profiles (temperature, humidity) from RAOB and a microwave
radiative transfer calculation. In step 4, the parameter α is calculated by using Eq. [4] for eachfrequency and each polarization. In step 5, the transmittance is calculated from the AMSR2
brightness temperature by using Eq. [5] for both polarizations. In step 6, the transmittances for eachfrequency are estimated as the geometric mean of transmittance found for each polarization. Then,
we repeat from step 3 to step 6 iteratively to find an optimized value for the atmospherictransmittance. The iterative algorithm is fairly stable, so that a few iterations are sufficient to reachthe convergence.
A further step of the algorithm consists in calculating two index parameters that can be converted
into TPW and CLW through a monotonic law. The indices, PWI (Precipitable Water Index) and CWI(Cloud liquid Water Index), are calculated by the following equation:

PWI = β ln(Tr19 HV ) − ln(Tr 22V ) + γ CWI[6]CWI = ln(Tr19 HV ) − ln(Tr 37 HV ) 2

[7]

2

2

2

where β = β (Ts ) is a function of SST and γ is a constant. The indices are defined to representmonotonically increasing (or decreasing) amounts with TPW and CLW increase (or decrease)(Takeuchi, 2002; Takeuchi et al., 2004). The third LUT, which is used in this conversion, defines therelationship between SST and the parameters β and γ . The LUT is produced in advance tomaximize correlations between PWI and observed RAOB TPW. An additional atmospheric moisturecorrection is applied to the CWI before it is being converted into CLW. As the CWI contains theatmospheric water vapor dependency in the intermediate parameters, the CWIs are stratified intosome TPW categories and the dependency is removed with a linear correction function. Pre-definedcoefficients are used in the conversion. The observation at 37 GHz contains some SSW dependency.
The SSW dependency is removed by using a wind speed related parameter (S36, defined in Shibata,2012) in the algorithm. The LUTs used in the algorithm are smoothed for possible input parameter

2-6

NDX-120015A

ranges to avoid discontinuities in the retrieval due to insufficient sampling data from RAOB.

4. Algorithm Implementation
4.1.

Implementation

The algorithm is implemented into the operational JAXA computer in the GCOM ground system todetermine TPW and CLW from GCOM-W1/AMSR2 data in real time. Additional quality flags
indicating the retrieval accuracy, as well as sea ice detection and rain flags, are generated. The
products are saved in files for every ascending and descending orbit. The data file format is HDF5.

4.2.

Input/output parameters

Input parameters
AMSR2 Level 1R brightness temperatures, i.e., measurements with dual polarization channels at 19
GHz, 23 GHz, 37 GHz, and 89 GHz, are used as input parameters for the algorithm. Information at
89 GHz is used mainly for rain detection by the algorithm. Several ancillary data (SST, SSW, andtemperature at 850 hPa) are interpolated spatially and temporally at the AMSR2 observation location.
The incidence angle of the AMSR2 measurement is used for the radiative transfer calculation in thealgorithm. As described in Section 2, three different look-up tables are used as input information by
the algorithm. The first one is for the computation of ocean emissivity. The inputs of the LUT are
SST, SSW, and AMSR2 channel microwave frequency (i.e., 19 GHz, 23 GHz, and 37 GHz). The
LUT produces ocean emissivity for each AMSR2 channel. The second one is for the computation of
atmospheric mean emission temperatures. The inputs of the LUT are the square of the atmospherictransmittance and the temperature at 850 hPa from ancillary data. The last one is for the conversion
from index parameters for atmospheric integrated water vapor and cloud liquid water to real values.
The LUTs are predefined by using a radiative transfer model and a collocated data set with AMSR2
and atmospheric profiles from RAOB.
Output parameters
The output parameters are atmospheric integrated water vapor and cloud liquid water in the
AMSR2 field of view. The algorithm produces some quality flags indicating the retrieval accuracy.
The data are categorized into three meteorological conditions (clear, cloudy, and rainy). Rain and sea
ice contamination flags are also produced.

4.3.

Ancillary data

The temperature at 850 hPa and SSW from JMA global atmospheric analysis, and SST from JMA
global sea surface temperature analyses, are used as ancillary data for the standard product. For realtime production of TPW and CLW, the analyzed values are not available. Temperatures at 850 hPa
and SSW from the short range global forecast of the JMA global model are used as ancillary input
2-7

NDX-120015A

data.

4.4.

Processing flow

Figure 1. Schematic diagram of the atmospheric water vapor (TPW) and cloud liquid water (CLW)
retrieval algorithm

2-8

NDX-120015A

4.5.

Example output

The global map of atmospheric integrated water vapor (Figure 2(a)) represents the moist region of
associated with the Inter-Tropical Convergence Zone (ITCZ) and the maritime continent. In the
global map of cloud liquid water (Figure 2(b)), cloudiness associated with low pressure systems isshown in the mid- and high latitude areas.
(a)

(b)

Figure 2. Global map of AMSR-E integrated water vapor (a) and cloud liquid water (b) on January 1,
2010 for ascending data. The units are [kg/m2] for water vapor and [kg/m2] for cloud liquid water.

4.6.

Limitations

Because the algorithm is applicable only to open ocean areas and rain free conditions, data overland and for rainy conditions are not available. Furthermore, as the algorithm does not include cloud
distribution inside the AMSR2 field of view, there is a possibility of overestimation orunderestimation of the cloud liquid water amount. It is difficult to validate the product withindependent real observations of cloud liquid water content.

5. Validation Concept
Verification of integrated water vapor
Historically, integrated water vapor calculated from radiosonde observation profiles has been usedas truth data in the verification of satellite microwave radiometer integrated water vapor products.
Since RAOB provides direct measurements of atmospheric humidity profiles, the estimatedintegrated water vapor can be considered highly accurate compared with other remote sensing
observations. In order to minimize land contamination in the integrated water vapor verification,
RAOB in small islands should be used in preference to stations on larger land masses, andappropriate temporal and spatial collocation criteria should be set to obtain sufficient sampling data
between GCOM-W1/AMSR2 and RAOB. In addition, careful quality control should be applied to
the radiosonde observations before the verification. However, the temporal and spatial resolution of
the RAOB network is not sufficient for the verification of satellite-based water vapor retrievals
under various meteorological conditions, due to the sparse observation network. In particular, the
2-9

NDX-120015A

local time of GCOM-W1/AMSR2 observations is fixed (13:30) as it is in the A-train orbit. This
limits the number of RAOBs available for the collocation process.
In addition to verification with RAOB, measurements from the Global Positioning System (GPS)have been made available to the public domain, and the Integrated Water Vapor (IWV) fromground-based GPS receivers can be used as new references for satellite product verification.
GPS IWV observations are also becoming valuable in operational NWP and weather forecasting,thanks to their high accuracy and frequent sampling (five-minute intervals). Preliminary validation
results for the AMSR-E TPW products are shown in Kazumori et al. (2012).

Verification of cloud liquid water
Since there are few direct CLW measurements available, it is difficult to verify the accuracy of the
CLW product comprehensively. Therefore, the accuracy of CLW products can be estimated only
indirectly. A proposed method in Greenwald et al. (2007) and Greenwald (2009) can be used for theestimation of certain CLW error characteristics. Cloud clear scenes of AMSR2 are selected using
MODIS collocated images. The cloud fraction of MODIS Level 2 products (Menzel et al., 2010) was
utilized in the cloud clear scene determination (i.e., zero cloud fraction as determined by collocated
MODIS products). The data are stratified with respect to SSW, SST, and TPW categories. This
technique has been used by Wentz and Meissner (2000) and enables us to identify false correlations
among different geophysical parameters derived from the microwave radiometer.

References
Greenwald, T. J., T. S. L’Ecuyer, and S. A. Chirstopher 2007: Evaluating specific error
characteristics of microwave-derived cloud liquid water products. Geophys. Res. Lett., 34,
L22807, doi:10.1029/2007GL031180.
Greenwald, T. J. 2009: A 2 year comparison of AMSR-E and MODIS cloud liquid water pathobservations. Geophys. Res. Lett., 36, L20805, doi:10.1029/2009GL040394.
Janssen, M A., Ed. (1993) Atmospheric Remote Sensing By Microwave Radiometry. John Wiley &Sons, Inc., pp. 572.
Japan Meteorological Agency 2007: Outline of Operational Numerical Weather Prediction at Japan
Meteorological Agency. Appendix to WMO technical progress report on the global dataprocessing and forecasting system and numerical weather prediction, Japan Meteorologial
Agency, Tokyo.
Available from http://www.jma.go.jp/jma/jma-eng/jma-center/nwp/outline-nwp/index.htm
Kazumori, M., Q. Liu, R. Treadon, and J. C. Derber 2008: Impact Study of AMSR-E Radiances in
the NCEP Global Data Assimilation System. Monthly Weather Review, 136, 541 – 559.

2-10

NDX-120015A

Kazumori, M., T. Egawa and K. Yoshimoto 2012: A retrieval algorithm of atmospheric water vapor
and cloud liquid water for AMSR-E. European Journal of Remote Sensing, 45, 63 – 74. doi:
10.5721/EuJRS20124507.
Menzel, M. P., R. A. Frey, and B. A. Baum 2010: Cloud top properties and cloud phase algorithm
theoretical

bases

document.

UniversityofWisconsin,

Madison.

Availablefrom

http://modis-atmos.gsfc.nasa.gov/_docs/CTP_ATBD_oct10.pdf.
Shibata, A. 2012: Ocean Wind Speed Retrieval Algorithm using the frequency 36GHz
Vertical/Horizontal and 6GHz Horizontal Data of the Advanced Microwave ScanningRadiometer (AMSR). European Journal of Remote Sensing, 45, 133 – 140. doi:
10.5721/EuJRS20124513.
Stephens G. L. 1994: Remote sensing of the lower atmospheric: an introduction. Oxford University
Press, pp. 523.
Takeuchi, Y. 2002: Algorithm theoretical basis document (ATBD) of the algorithm to derive total
water vapor content from ADEOS-II/AMSR. EORC Bulletin/Technical Report.
Takeuchi, Y., T. Tauchi, S. Saito, K. Imaoka, M. Kachi, and A. Shibata 2004: A total column
precipitable water algorithm for ADEOS-II/AMSR and Aqua/AMSR-E. Italian Journal of
Remote Sensing 30/31, 143 –157.
Wentz, F. J., and T. Meissner 2000: AMSR Ocean Algorithm Theoretical Basis Document (ATBD).
RSS Tech. Doc. 121599A –1, Remote Sens. Syst., Santa Rosa, Calif.

2-11

NDX-120015A

Chapter 3.
Description of GCOM-W1 AMSR2Precipitation Algorithm

Kazumasa Aonashi
Meteorological Research InstituteJapan Meteorological Agency
3-1

NDX-120015A

Table of Contents

1.

Introduction ................................................................................................................... 3

2.

Algorithm Overview ....................................................................................................... 4

3.

Theoretical Description .................................................................................................. 5

4.

Algorithm Implementation ............................................................................................. 9

5.

4.1.

Implementation ......................................................................................................... 9

4.2.

Input/output parameters ............................................................................................. 9

4.3.

Ancillary data ........................................................................................................... 9

4.4.

Processing flow ...................................................................................................... 10

4.5.

Example output....................................................................................................... 13

4.6.

Limitations ............................................................................................................. 18

Validation Concept ....................................................................................................... 18

References ........................................................................................................................... 18

3-2

NDX-120015A

1. Introduction
Satellite passive microwave radiometer (MWR) brightness temperatures (TBs)
give information on emission from raindrops and scattering by frozen particles.
MWRs recently have become the principal sensors for global precipitation
retrieval, since these emission and scattering signals (increased
lower-frequency TB and decreased higher-frequency TB, respectively, from a
precipitation-free atmosphere) have a more direct relationship with precipitation
rates than infrared radiometer (IR) cloud-top TBs (Ebert and Manton 1998).
MWR precipitation retrieval algorithms are generally classified as either
statistical or physical. Statistical algorithms (e.g., Grody 1991; McCollum andFerraro 2003) require observed precipitation data sets for derivation of theempirical relationship between precipitation rates and TBs. Physical algorithms(e.g., Panegrossi et al. 1998; Kummerow et al. 2001; Aonashi and Liu 2000)
require a priori models for precipitation-related variables, and radiative-transfer
models (RTMs) in order to calculate look-up tables (LUTs) relating toprecipitation rates and TBs.
Before the Tropical Rainfall Measuring Mission (TRMM) launch, global
information on precipitation-related variables was sparse. Accordingly,
conventional statistical algorithms depended on observed precipitation data setstaken over limited areas. Conventional physical algorithms resorted toprecipitation-related variable models based on a small number of field
campaigns (Aonashi and Liu 2000, hereafter referred to as AL) or
cloud-resolving model outputs (Panegrossi et al. 1998; Kummerow et al. 2001).
Recent studies using TRMM data, however, have revealed the global
distribution of various precipitation-related variables, such as particle size andprecipitation types (e.g., Kozu et al. 2009; Takayabu 2006). The TRMM data
have also enabled us to research the statistical relationship between TBs andprecipitation-related variables (e.g., Seto et al. 2005; Kubota et al. 2007, 2009;
Kida et al. 2009).
Under the Global Satellite Mapping of Precipitation project (GSMaP), we
developed a new physical precipitation retrieval algorithm (hereafter referred toas the GSMaP algorithm) by improving the AL algorithm based on the above
studies.
The AL algorithm finds the surface precipitation rates that give the best fit for
the forward-calculated drop in TB at 37 GHz (TB37) and 85 GHz (TB85) with the

3-3

NDX-120015A

MWR observations. The AL algorithm underestimates the TRMM Precipitation
Radar (PR) surface rain (Rainsurf), especially for shallow precipitation events.
This algorithm also shows differences in retrieval biases between stratiform andconvective precipitation.
We introduced indices for the frozen precipitation depth and stratiform rain
ratio from MWR TBs into the retrieval algorithm in order to reduce the above
systematic over-land errors.
The purpose of the present study is to describe the GSMaP algorithm and toreport its validation results. To this end, we compared its retrievals from TRMM
Microwave Imager (TMI) TBs in 1998 with the TRMM PR data.

2. Algorithm Overview
The underlying concept of the precipitation retrieval algorithm is to findprecipitation rates that give Radiative Transfer Model (RTM)-calculated TBs thatbest fit the MWR TBs. This algorithm employs Polarization Corrected
Temperatures (PCT) at higher frequencies (36 and 89 GHz for AMSR2, 37 and
85 GHz for TMI, hereafter referred to as 37 and 85 GHz, for simplicity) over landand coast, and TBs with vertical polarization at lower frequencies (10, 19, and 36
GHz for AMSR2; 10, 19, and 37 GHz for TMI, hereafter referred to as 10, 19, and
37 GHz) in addition to PCT at 37 GHz (PCT37) and 85 GHz (PCT85) over ocean.
The precipitation retrieval algorithm consists of the forward calculation part tocalculate the Look-Up Tables (LUTs), and the retrieval part to estimate
precipitation rates from the observed TBs using the LUTs.
In the forward calculation part of the algorithm, we adopted atmospherictemperature, Freezing Level Height (FLH), and surface temperature over landand coast from JMA Global Analysis data. We also employed precipitation
profiles and raindrop size distribution models based on the TRMM observation
studies. We calculated LUTs by incorporating the above atmospheric andprecipitation variables into the RTM.
In the retrieval part of the algorithm, we first detected precipitation areas andestimated precipitation inhomogeneity based on the TRMM observation studies.
Then we retrieved precipitation using the scattering signals. As scattering
signals, we choose PCT37 in addition to PCT85, in spite of the coarser
resolution and smaller signals of PCT37. This decision is based on a statisticalcomparison between the TRMM Microwave Imager (TMI) scattering signals and3-4

NDX-120015A

TRMM Precipitation Radar (PR) precipitation. For retrieval over ocean, emission
signals at lower frequencies were available, in addition to the scattering signals.
The retrieval algorithm found the optimal precipitation over ocean by minimizinga cost function that was a weighted sum of squares of differences between theobservation and the forward calculation for TB10v and TB19v.

3. Theoretical Description
3.1 Precipitation-related variables
The GSMaP algorithm consists of a forward-calculation part to calculate the
LUTs and a retrieval part to estimate precipitation rates from the observed TBs
using the LUTs.
As Figure 1 schematically illustrates, the RTM calculation requires information
on atmospheric variables (temperature, relative humidity, freezing level height(FLH), cloud liquid water content (CLWC), etc.), surface variables (surface
temperature, sea-surface wind speed), and precipitation-related variables(precipitation profiles, precipitating-particle size distribution, refractivity of frozen
and mixed-phased particles, horizontal inhomogeneity of precipitation, etc.).
This subsection will describe the above physical variables and the RTM
calculation methods used in the GSMaP algorithm. As described in the
Introduction, the main improvement in the forward calculation part over AL is theuse of new precipitation-related variable models, as summarized in Table 1.
Fig. 1. Schematic illustration of information required for the RTM calculation.

3-5

NDX-120015A

Table 1. Main differences in precipitation-related variable models of the AL and GSMaP
algorithms
Precipitation-related variables
Precipitation profile
Raindrop size distribution
Mixed-phase particle size
distribution and refractivity

AL
Constant vertical gradients
Marshall-Palmer distribution
Neglected

GSMaP
The PR profiles averaged foreach precipitation typeKozu et al. (2008) model
Nishitsuji et al. (1983) model
for stratiform precipitation

We adopted the following convective and stratiform precipitation models forthe precipitation-related variables:(1) Precipitation Profile
We constructed statistical precipitation profile models using TRMM
observations. To this end, we first classified PR data into 10 types (six over land,
four over ocean and coast) using the PR precipitation parameters (rain area,
stratiform rain-area fraction, precipitation-top level, etc.) and the ratio between
the PR precipitation rates and TRMM Lightning Imaging Sensor (LIS) flash rates(Takayabu 2006).
We then produced convective and stratiform precipitation profile models forthese types by averaging the PR convective and stratiform precipitation profiles
over prescribed precipitation ranges for each type. In this averaging, we used
3-6

NDX-120015A

profiles relative to FLH in order to exclude the influence of atmospherictemperature variations (Kubota et al. 2007).
(2) Precipitating-particle size distribution models
For the raindrop-size distribution, we constructed statistical models for eachprecipitation type using globally available PR observations (Kozu et al. 2009). To
this end, we first assumed a gamma distribution for the raindrop size:=
N ( D) N 0 D µ exp(−λ D) ,（１）

where N is the number concentration for particles with diameter D, µ = 3, andN 0 and λ are parameters to be determined. For convective precipitation, we
calibrated N 0 and λ using attenuation-corrected factors derived from the PR
data. For stratiform precipitation, we used the standard values for N 0 and λ of
PR 2A25 version 6. We then averaged the above values for each precipitation
type to produce statistical models of raindrop-size distribution.
Conventional models were used for frozen and mixed-phase particle-size
distributions that could not be estimated from the TRMM observations. For thefrozen-particle size distribution, we applied the Marshall–Palmer distribution toboth convective and stratiform precipitation, similar to AL. The particle-size
distribution for mixed-phase stratiform precipitation (between FLH minus 1 km
and FLH) was parameterized in terms of atmospheric temperature (Nishitsuji etal. 1983). We neglected mixed-phase convective precipitation.
(3) Refractivity of frozen and mixed-phase particles
We calculated the refractivity of convective and stratiform frozen particles,
assuming them to be a mixture of ice and air with an empirically prescribed
constant density (200 kg m-3), similar to AL. The refractivity of mixed-phasestratiform particles was parameterized in terms of atmospheric temperature(Nishitsuji et al. 1983).
3.2 Indices of the frozen precipitation depth and stratiform rain ratio
from MWR TBs
We introduced the ratio of decreases in TB85 to decreases in TB37 (R8537)
as the index of the frozen precipitation depth. (We expressed R8537 in terms ofratio of precipitation retrieved from the decrease in TB85 (Rain85) to those from
3-7

NDX-120015A

TB37 (Rain37) using the conventional GSMaP algorithm.) As the index ofstratiform rain ratio, we introduced the horizontal precipitation inhomogeneity
derived from Rain85 (Sigma85). Then we classified the TMI retrievals with
R8537 and Sigma85, and compared them with PR Rainsurf for 1998. The resultsshow:1) Both Rain37 and Rain85 greatly underestimated Rainsurf for cases with
R8537 <0.7 (See Fig. 2 and Fig. 3).
2) Precipitation with Sigma85 <1 tended to have larger scattering for caseswith R8537 >0.7, compared to those with Sigma85 >1.
3) The above relations had little seasonal variation and held good for all
precipitation types, except for high land and wintertime extra-tropical low
types that showed a large underestimation of Rainsurf for all R8537 andSigma85 ranges.
Fig. 2. Scatter diagrams between TRMM PR Rainsurf and TMI Rain37 overland for 1998, classified with R8537 and Sigma85.

Fig. 3. Scatter diagrams between TRMM PR Rainsurf and TMI Rain85 overland for 1998, classified with R8537 and Sigma85.

3-8

NDX-120015A

4. Algorithm Implementation
4.1. Implementation

The underlying concept of the GSMaP algorithm is to determine theprecipitation rates that yield RTM-calculated TBs that best fit with MWR TBs.
The GSMaP algorithm employs PCT37 and PCT85 over land and coast and TBs
with vertical polarization at 10, 19, and 37 GHz (TB10v, TB19V, and TB37v), in
addition to PCT37 and PCT85 over ocean. The GSMaP algorithm consists of aforward-calculation part to calculate the LUTs and a retrieval part to estimate
precipitation rates from the observed TBs using the LUTs.
4.2. Input/output parameters

The algorithm uses AMSR2 TB as its input.
4.3. Ancillary data

We adopted the atmospheric temperature, FLH, and surface temperature overland and coast from the Japan Meteorological Agency (JMA) global analysis(GANAL) data, since GANAL gives a priori, global information on these variables
with certain accuracy. Similarly, we adopted the sea-surface temperature fromthe JMA merged satellite and in situ data on global daily sea-surface

3-9

NDX-120015A

temperatures in the ocean (MGDSST). For relative humidity and CLWC, we
assumed the same constant values, 100% and 0.5 kg m-2, as AL.
The LUTs were calculated for every 5 × 5 degree point within thelatitude-longitude coordinates, as we will describe in Subsection 3.1.c. Thus, weincorporated the above variables for the 5 × 5 degree points into the RTM
calculation.
4.4. Processing flow
（１）Forward calculation
First, we calculated the convective and stratiform LUTs for horizontallyhomogeneous precipitation by incorporating the above atmospheric-, surface-,and precipitation-related variables into the RTM program of Liu (1998). This
program computed TBs for a plane-parallel atmosphere using a four-stream
approximation. In this program, all precipitation particles were assumed to be
spherical, and the absorption and scattering coefficients and the phase functions
were computed based on the Mie theory.
We then derived convective and stratiform LUTs with various horizontal
precipitation inhomogeneities using the approximation of AL. This approximationhconverted TBs for homogeneous precipitation ( TB ) into those foriinhomogeneous precipitation ( TB ), assuming a lognormal distribution of the
precipitation:11iexp[− 2 (ln (pr ) − ln (pr )) 2 ]TB h d ln (pr ),TB = ∫
2x
2p x

(2)

where ln(pr) is the natural logarithm of the precipitation, and ln (pr ) and ξ
are the average and the standard deviation of ln(pr).
The LUTs used in the retrieval part were weighted averages of the above
convective and stratiform LUTs. The weights were determined by the statisticalfrequency distribution of the PR convective and stratiform precipitation for eachprecipitation type and surface precipitation rate.
（２）Retrieval part
This subsection will describe the procedures followed in the retrieval part of
the algorithm, schematically illustrated in Fig. 4. As described in the Introduction,
the main improvements in the retrieval part over AL are as follows: (1) new
precipitation detection and inhomogeneity estimation methods based on TRMM
3-10

NDX-120015A

observation studies; and (2) a new scattering retrieval method using PCT37 and
PCT85 and scattering-signal correction for tall precipitation.
Fig. 4. Schematic illustration of procedures in the retrieval part of the GSMaP
algorithm.

(a) Detection of precipitating areas
Over land, we adopted the detection methods of Seto et al. (2005). We
calculated the scattering indices from TMI TB85v and TB21v, and judged thatprecipitating areas had scattering indices larger than the statistically determinedthresholds. Since this method mistook snow cover for precipitation, we regarded
all areas with surface temperatures lower than 273.2 K as precipitation-free.
Over coastal areas, we used the method of Kubota et al. (2007). This is an
improvement of McCollum and Ferraro (2005) that detects precipitating areasusing PCT85 scattering signals and a decision tree of several empiricalconditions for TMI TBs.
Over ocean, we adopted the method of Kida et al. (2009). We judged that
3-11

NDX-120015A

precipitating areas had either PCT85 scattering or TB37v emission signals,
which were expressed by the respective decreases in TMI PCT85 and increases
in TB37v from the LUT TBs for 0 mm hr-1.
(b) Estimation of precipitation inhomogeneity
We used high-resolution PCT85 to estimate the standard deviation of thenatural logarithm of the precipitation (SDNLPR), similar to AL.
First, we converted the PCT85 scattering signals into precipitation (Rain85)
using LUTs for a prescribed value of SDNLPR (=1.0), and calculated a first
guess of the inhomogeneity (Sigma85) from Rain85 within the TB10v field ofviews (FOVs).
This first guess was then adjusted based on the statistical relationshipbetween Sigma85 and SDNLPR calculated from PR (Kubota et al. 2009). The
LUTs for this adjusted value were used in the following retrieval procedures.
(c) Precipitation retrieval using scattering signals
As scattering signals for precipitation retrieval, we chose PCT37 in addition to
PCT85, in spite of the coarser resolution and smaller signals of PCT37.
We introduced R8537 as the index of the frozen precipitation depth. As theindex of stratiform rain ratio, we introduced Sigma85. Then, we classified the
TMI retrievals with R8537 and Sigma85, and compared them with PR Rainsurf
for 1998. Next, we derived linear fitting coefficients between Rain37, Rain85,
and Rainsurf for each R8537 and Sigma85 class and precipitation type for 1998.
The retrieval part used these fitting coefficients for the calibration.
(d) Precipitation retrieval using emission signals over ocean
For retrieval over ocean, emission signals at lower frequencies were available
in addition to the scattering signals. Similar to AL, the GSMaP algorithm foundthe optimal precipitation over ocean by minimizing a cost function (J) that was aweighted sum of squares of the differences between the observation and theforward calculation for TB10v and TB19v:
km
2
1
(3)
J ∑∑ 2 TB of ( yˆ k ) − TB if  yˆ k : r ( xˆ j )  ,f k =1 2σ f{}where r ( xˆ j ) is the precipitation at point xˆ j , yˆ k is the location of a TMI pixel,

3-12

NDX-120015A

and TB of and TB if

represent the TMI TBs and the FOV-averaged TBs

calculated from r ( xˆ j ) at frequency f. σ f is the standard deviation of thedifferences between TB of and TB if . Here, ∑ represents summation for f = 10f

and 19 GHz, and ∑

km
k =1

is the summation of TMI pixels over ocean.

We derived the first guess for this minimization by adjusting the scattering
retrievals based on precipitation retrieved from TB10v (Rain10v). In thisadjustment, we multiplied the scattering retrievals by the 10 GHz FOV averagesof the ratio of scattering retrievals to Rain10v. For shallow precipitation with zeroprecipitation retrievals, we employed precipitation retrieved from TB37v(Rain37v) as the first guess.
We obtained the optimal precipitation over ocean by AL’s method of solving
the gradient equation of the cost function.
4.5 Example output

We used the GSMaP algorithm for retrieving TMI global precipitation
distributions for January, April, July, and October 1998. We then validated theseprecipitation retrievals using the PR precipitation rates.
Fig. 5. Global precipitation distribution of the GSMaP TMI retrieval for (a)January, (b) April, (c) July, and (d) October 1998.
(a)

3-13

NDX-120015A

(b)

(c)

(d)

3-14

NDX-120015A

Fig. 6. Global precipitation distribution of the TRMM PR Rainsurf for (a)January, (b) April, (c ) July, and (d) October 1998.
(a)

(b)

(c)

3-15

NDX-120015A(d)Fig. 7. Scatter diagrams between GSMaP TMI retrieval and PR Rainsurf oversea for (a) January, (b) April, (c) July, and (d) October 1998.
(a)

(b)

3-16

NDX-120015A

(c)

(d)

Fig. 8. Scatter diagrams between GSMaP TMI retrieval and PR Rainsurf overland for (a) January, (b) April, (c) July, and (d) October 1998.
(a)
(b)

(c)

(d)

3-17

NDX-120015A

4.6 Limitations

5. Validation Concept(TBD)

References
Aonashi, K. and G. Liu, 2000: Passive Microwave Precipitation Retrievals using TMI during
the Baiu period of 1998. Part 1: Algorithm description and Validation. J. Appl. Meteor. 39,
2024-2037.
Ebert, E.E., and M.J. Manton, 1998: Performance of satellite rainfall estimation algorithms
during TOGA COARE. J. Atmos. Sci., 55, 1537-1557.
Grody, N., 1991: Classification of snow cover and precipitation using the Special Sensor
Microwave Imager, J. Geophys. Res., 96(D4), 7423-7435.
Iguchi, T., 2007: Space-borne radar algorithms in measuring precipitation fromspace—EURAINSAT and the future, edited by V. Levizzani, P. Bauer, and F.J. Turk,
Springer, 199-212.
Iguchi, T., T. Kozu, J. Kwiatkowski, R. Meneghini, J. Awaka, and K. Okamoto, 2009:

Uncertainties in the rain profiling algorithm for the TRMM precipitation radar. J. Meteor.
Soc. Japan, this issue.
Kida, S., S. Shige, T. Kubota, K. Aonashi, and K. Okamoto, 2009: Improvement ofrain/no-rain classification methods for microwave radiometer observations over oceanusing the 37-GHz emission signature. J. Meteor. Soc. Japan, this issue.
Kozu, T., T. Iguchi, T. Kubota, N. Yoshida, S. Seto, J. Kwiatkowski, and Y. N. Takayabu,
2009: Feasibility of raindrop size distribution parameter estimation with TRMM
precipitation radar. J. Meteor. Soc. Japan, this issue.
Kubota, T., S. Shige, H. Hashizume, K. Aonashi, N. Takahashi, S. Seto, M. Hirose, Y.N.
Takayabu, K. Nakagawa, K. Iwanami, T. Ushio, M. Kachi, and K. Okamoto, 2007: Global
precipitation map using satellite-borne microwave radiometers by the GSMaP project :production and validation. Trans. Geosci. Remote Sens., 45, 7, 2259-2275.
Kubota, T., S. Shige, K. Aonashi, and K. Okamoto, 2009: Development of nonuniform
beam-filling correction method in rainfall retrievals for passive microwave radiometers
over ocean using TRMM observations. J. Meteor. Soc. Japan, this issue.
Kummerow, C., Y. Hong, W.S. Olson, S. Yang, R.F. Adler, J. McCollum, R. Ferraro, G. Petty,
D.B. Shin, and T.T. Wilheit, 2001: The Evolution of the Goddard profiling algorithm
3-18

NDX-120015A(GPROF) for rainfall estimation from passive microwave sensors. J. Appl. Meteor., 40,
1801-1820.
Kummerow, C., P. Poyner, W. Berg, and J. Thomas-Stahle, 2004: The effects of rainfall
inhomogeneity on climate variability of rainfall estimated from passive microwave
sensors. J. Atmos. Oceanic Technol., 21, 624-638.
Liu, G., 1998: A fast and accurate model for microwave radiance calculation. J. Meteor. Soc.
Japan, 76, 335-343.
McCollum J.R., and R.R. Ferraro, 2003: Next generation of NOAA/NESDIS TMI, SSM/I, and
AMSR-E microwave land rainfall algorithms. J. Geophys. Res., 108 (D8), 8382,
doi:10.1029/2001JD001512.
McCollum, J.R., and R.R. Ferraro, 2005: Microwave rainfall estimation along coasts. J.
Atmos. Oceanic Technol., 22, 497-512.
Nishitsuji, A., M. Hoshiyama, J. Awaka, and Y. Furuhama,1983: An analysis of propagative
character at 34.5 GHz and 11.5 GHz between ETS-II satellite and Kasima station –On the
precipitation model from stratus. IEICE Trans. (in Japanese), J66-B, 1163-1170.
Olson, W.S., D.K. Christian, S. Yang, G.W. Petty, W.-K. Tao, T.L. Bell, S.A. Braun, Y. Wang,
S.E. Lang, D.E. Johnson, and C. Chiu, 2006: Precipitation and latent heating distributions
from satellite passive microwave radiometry. Part I: Improved method and uncertainties. J.
Appl. Meteor. Climat., 45, 702-720.
Panegrossi, G., S. Dietrich, F.S. Marzano, A. Mugnai, E.A. Smith, X. Xiang, G.J. Tripoli, P.K.
Wang, and J.P.V. Poiares Baptista, 1998: Use of cloud model microphysics for passive
microwave-based precipitation retrieval: significance of consistency between model andmeasurement manifolds. J. Atmos. Sci., 55, 1644-1673.
Seto, S., N. Takahashi, and T. Iguchi, 2005: Rain/no-rain classification methods for
microwave radiometer observations over land using statistical information for brightness
temperatures under no-rain conditions. J. Appl. Meteor. 44. 1243-1259.
Shige, S., H. Sasaki, K. Okamoto, and T. Iguchi, 2006: Validation of rainfall estimates fromthe TRMM precipitation radar and microwave imager using a radiative transfer model: 1.
Comparison of the version-5 and -6 products. Geophys. Res. Lett., 33, L13803,
doi:10.1029/2006GL026350.
Shige, S., T. Watanabe, H. Sasaki, T. Kubota, S. Kida, and K. Okamoto, 2008: Validation ofwestern and eastern Pacific rainfall estimates from the TRMM PR using a radiative transfer
model. J. Geophys. Res., 113, D15116, doi:10.1029/2007JD009002.
Takayabu N. Y., Rain-yield per flash calculated from TRMM PR and LIS data and its
relationship to the contribution of tall convective rain. Geophys. Res. Lett., 33, L18705,
doi:10.1029/2006GL027531, Sept 2006.
3-19

NDX-120015A

Chapter 4.
Description of GCOM-W1 AMSR2
Sea Surface Temperature Algorithm

Akira Shibata
Meteorological Research InstituteJapan Meteorological Agency

4-1

NDX-120015A

Table of Contents
1.

Introduction

3

2.

Algorithm Overview

3

3.

Theoretical Description

4

4.

Algorithm Implementation

105.

4.1 Implementation

10

4.2 Input/output parameters

10

4.3 Ancillary data

10

4.4 Processing flow

10

4.5 Example output

11

4.6 Limitations

11

Validation Concept

11

References

4-2

NDX-120015A

1. Introduction
Passive microwave radiometers are capable of measuring sea surface temperature (SST) undercloud cover. Usable frequencies for measuring SST range from 4 to 10 GHz. It has been noticed thatthe sensitivity to low SST becomes less as the frequency is increased (Hollinger and Lo, 1984). The
first spaceborne microwave radiometer capable of measuring SST was the Scanning Multichannel
Microwave Radiometer (SMMR) carried on Seasat 1 and Nimbus 7, both of which were launched in
1978. Although SST retrievals from the SMMR on Nimbus-7 had an error as high as 1.12°C, most of
the SST error was due to calibration problems unique to the SMMR design (Milman and Wilheit,1985). The second radiometer capable of measuring SST was the TRMM Microwave Imager (TMI)
on the Tropical Rainfall Measuring Mission (TRMM) launched in 1997. The SST accuracy retrieved
from TMI at 10 GHz is much better than that of the SMMR (Shibata et al., 1999), and the citedresults seem to promise accurate SST retrievals at 6 GHz with future radiometers. The third sensors
were two AMSRs; one is the AMSR-E onboard the NASA Aqua launched in May 2002, and theother is the AMSR onboard the ADEOS-II launched in December 2002. The AMSRs were the first
sensors to provide an accurate microwave SST for the global ocean, and the averaged rms of thedifference between AMSR and buoy SSTs was 0.59°C for AMSR-E over an 18 month period, and
0.74°C for AMSR over seven months (Shibata, 2004). This paper will describe the SST algorithm
for AMSR2, onboard the GCOM-W1 satellite, which is the first satellite of Japan’s Global Change
Observation Mission.

2. Algorithm Overview
We retrieve the SST using the brightness temperatures at 6 GHz with vertical polarization(V-pol), referred to henceforward as 6V. In addition to the SST, 6V contains various signals related
to other parameters. These are (a) atmospheric effect, (b) wind effect, (c) salinity effect, (d) land
contamination, (e) sea ice contamination, and (f) sun glitter contamination. Among the first threeeffects, (c) is very small. We can correct (c) by using the monthly climate salinity. Effects (a) and (b)are large, and we will explain in detail methods for correcting them. Contaminations due to (d) to (f)are also significant, so we eliminate contaminated areas as much as possible. In our algorithm, areas
missing owing to (d) are within about 100 km of the shorelines, and areas missing owing to (f) occur
where there is an angle of 25° or less between the AMSR viewing direction and the sun glitter
direction. In (e), the sea ice case, we eliminate those pixels where sea ice is detected by AMSR.
Effects (a) to (c) and contaminations (d) to (f) are common to microwave radiometers in general. In
addition to these, we need to correct other errors specific to AMSR. These are (g) incidence angle

4-3

NDX-120015A

variation, (h) scan bias error, and (i) radio frequency interference (RFI). In the case of AMSR, themain cause of (h) is a cold sky mirror blocking a main reflector viewing at the rightmost scanning
position. A temperature drop of −1.5 K may be registered over the ocean by 6V when scanning at
the rightmost position (for a 196 pixels swath). At frequencies other than 6 GHz, no scan bias isdiscernible. RFI at 6.9 GHz is very severe over continents and sometimes around islands in theocean.
In our algorithm to retrieve the SST, we require an initial SST value as one of the inputs. We
specify this initial value from the Reynolds SST for AMSR and AMSR-E, and from the Merged
satellite and in-situ data Global Daily SST (MGDSST) of the Japan Meteorological Agency (JMA)for AMSR2.

3. Theoretical Description
Atmospheric Corrections
The atmospheric effects at 6 GHz are mainly attributed to molecules (water vapor and oxygen)and cloud liquid water. These effects can be estimated by using 23V and 36V (V-pol of brightness
temperature at 23 GHz and 36 GHz, respectively). The 23 GHz channel is sensitive to water vapor(there is an absorption line of water vapor at 22.235 GHz), while the 36 GHz channel is sensitive tocloud liquid water. Only V-pol is used because it is less sensitive to ocean wind than the horizontal
polarization (H-pol) (Shibata, 2003). Fig. 1 shows an example of atmospheric correction for 6V(units: K) at an SST of 20°C, where the horizontal axis is 36V and the vertical axis is 23V. To makethe atmospheric corrections, we subtract these values from the AMSR 6V.

Fig. 1 Atmospheric correction value for 6V

4-4

NDX-120015A

We obtained Fig. 1 by using a microwave transfer model typically applied to one-year Japanese
aerological observations (Shibata, 1994). The microwave model consists of 60 levels ranging fromthe surface to 12 km, each 200 m in depth, followed by 18 levels for each 1 km of depth from 12 to
30 km. Within each level, the microwave signal decays owing to absorption by molecules and cloudliquid water. Molecules and cloud liquid water emit microwaves corresponding to the ambient
temperature. In the first step of the sequential calculations, an external microwave of 2.7 K arrives
from the cosmos, and the down-welling microwave is calculated for each level in the manner
described above. At the ocean surface, the microwave is reflected, with a reflection intensity given
by Fresnel’s formula (r). The ocean surface emits microwave radiation equal to SST × (1 - r). For 6V,
the value of r depends on SST and is roughly 0.6. The microwave signal in an upward direction fromthe ocean surface is calculated in the same manner as for the downward direction. Finally, themicrowave radiation at 30 km is obtained. The values shown in Fig. 1 were obtained by subtracting
the ocean surface radiation (i.e., SST × (1 - r)) from the values calculated at a height of 30 km; this isexpressed by Eq. (1).
atmos_effect_6V(H) = simu_6V(H) - calm_ocean_6V(H),(1)

where calm_ocean_6V(H) is given by the Fresnel formula.
Japanese aerological stations are located between the sub-tropical and sub-arctic regions, and theaccumulated water vapor content calculated from this dataset ranges from 3 to 60 kg/m2 (or mm).
Values for 23V vary from 190 K to 260 K owing to these variations in water vapor. The cloud liquidwater is inserted into the microwave transfer model at heights between 1 and 5 km, where its values
range between 0 and 2 kg/m2. Values for 36V vary from 205 K to 260 K, mainly because of

4-5

NDX-120015A

variations in cloud liquid water.
Atmospheric effects on 6V are calculated from 0°C to 35°C SST at intervals of 5°C, and areinterpolated for the intermediate SSTs. Fig. 1 applies to 6V, and we have a similar figure for 6H. In
rainy conditions, the microwaves at 36 and 23 GHz are scattered by raindrops, and we are unable touse those data quantitatively without considering these scattering effects. In our algorithm, the SST
is set as missing under rainy conditions, defined as occurring when the atmospheric effects on 6V
are larger than 6.6 K. (We determined the value of 6.6 K from the actual AMSR data.)

Wind Effect Corrections
Wind effects have two components: wind speed and wind direction. To observe these effects at
6 GHz, we created Fig. 2 by combining data from AMSR and SeaWinds on ADEOS-II. Fig. 2 (a)shows the relation of 6V to relative wind direction (RWD); (b) shows the relation of 6H to RWD,
and (c) shows the relation between 6V and 6H. The units of the vertical axis in Figs.2 (a) and (b),and of both axes in Fig. 2 (c), are Kelvin, annotated as 6V (K). To create Fig. 2, we used two days of
AMSR and SeaWinds data for the global ocean, for May 15 and 19, 2003. The RWD is defined as
the relative angle between AMSR viewing direction and the SeaWinds wind direction. An RWD of
180° corresponds to an upwind direction; 90° and 270° correspond to a crosswind direction; and 0°(and 360°) correspond to a downwind direction. SeaWinds wind speeds from 1 to 17 m/s at 4 m/s
intervals are shown in Fig. 2.
6V and 6H in Fig. 2 were already corrected for the atmospheric effect and SST effect, as
expressed in Eq. (2).
6V(H) = AMSR_6V(H) – atmos_effect_6V(H) – calm_ocean_6V(H),(2)

Fig. 2 Characteristics of (a) 6V and (b) 6H as a function of wind speed and relative direction(a) 6V versus RWD

4-6

NDX-120015A(b) 6H versus RWD(c) 6V versus 6H

4-7

NDX-120015A

In Fig. 2 (a), it can be seen that 6V does not change with RWD for wind speeds of 1 and 5 m/s,
but that it increases at higher wind speeds. 6V varies with RWD at wind speeds exceeding 9 m/s,
reaching a maximum in the upwind direction and a minimum in the downwind direction. The
difference between the maximum and minimum values is 1.5K at a wind speed of 13 m/s. In Fig. 2(b), it can be seen that 6H increases even at small wind speeds, such as 1 m/s. 6H varies with therelative wind direction at wind speeds higher than 9 m/s, reaching a maximum in the crosswind
direction and a minimum in the upwind and downwind directions.
In Fig. 2 (c), it can be seen that 6V does not change over a range where 6H is less than 3.8K =z0, and that both 6V and 6H increase above z0. The slope of 6V/6H (designated as sp) ranges from
0.5 to 0.7 (0.5 corresponds to the downwind direction; 0.7 corresponds to the upwind direction;
values for the crosswind direction are intermediate). From Fig. 2 (c), we can calculate the extent of
the wind effect on 6V (designated as inc_6V) as shown in Eq. (3).
inc_6V=0for 6H less than z0= (6H – z0) × sp

for 6H greater than z0(3)

To avoid SST errors under very high wind speeds, we regard the SST as missing when 6H
exceeds (z0 + 9).
Fig. 2 (c) is obtained from data for the global ocean. The relation between 6V and 6H depends
on the air-sea temperature difference (ASD), in addition to the RWD. To observe the dependency on
the ASD, we use air temperature data from the Global Analysis (GANAL) weather forecast model
operated by the Japan Meteorological Agency. Fig. 3 shows the 6V-6H relation in the downwind
direction and for three ASD cases (Shibata, 2007). The solid line corresponds to the global oceancase, the dotted line labeled by sp1 corresponds to the unstable case (i.e., ASD is –3°C), and thedotted line labeled by sp2 corresponds to the stable case (i.e., ASD is 2°C). The 6V-6H relation with
the ASD in the crosswind direction is similar to the one in the downwind direction; however, the
6V-6H relation with the ASD in the upwind direction is different. The slope is almost the same forthe three cases.

4-8

NDX-120015A

Fig. 3 Different slopes for different cases of ASD

To determine sp, we must know the RWD. We can determine the RWD from the AMSR data.
Fig. 4 shows the relation between 6H and S36. The parameter S36 represents the wind effect on the
36 GHz channel, and is described in detail in the SSW algorithm document.

Fig. 4 Method for Determination of RWD

In Fig. 4, it is seen that the data points for upwind, downwind, and crosswind RWD become
separated from one another for 6H above around 4 K. Below 4 K, the data points for the three
directions overlap one another.
The results from Fig. 4 indicate that we can determine the RWD by combining 6H and S36, i.e.,
we measure the distance from a point of (6H, S36) to the fitting curve corresponding to thecrosswind direction, and normalize it to a value between −1 and 1 (designated as dd). A value of −1
corresponds to upwind, 0 to crosswind, and 1 to downwind. The coefficient sp ranges between 0.5

4-9

NDX-120015A

and 0.7. Thus, it can be expressed as shown in Eq. (4),sp = 0.57 – 0.13 × dd
= 0.57 – 0.07 × dd

−1 < dd ≦0

for upwind

0 ≦ dd < 1.0 for downwind(4).

Conversion to SST
At this point, we have already corrected 6V for atmospheric and wind effects. We now retrieve
the SST from the corrected 6V by using the relation between 6V and SST. Fig. 5 shows this relation.
The curves shown in Fig. 5 were obtained by using the Fresnel formula, with an incidence angle of
55.0° and a salinity of 3.5%. The complex dielectric constant of the ocean water was obtained from
Klein and Swift (1977). Since the complex dielectric constant seems not to be quite perfect, we madesome small adjustments to the theoretical values for 6V. These adjustments became larger in caseswith lower SST values (e.g., −0.5 K at SST = 0°C; 0.0 K at SST = 30°C).
In most cases, the AMSR SST is obtained from 6V. In rare cases, such as in areas where RFI isexperienced at 6 GHz, the AMSR SST was obtained from the vertical polarization of the brightness
temperature for the 10 GHz channel (10V). The 10V SST was retrieved only at 10°C and higher,
because the sensitivity of 10V to SST is very weak below 10°C, as seen in Fig. 5. The atmospheric
and wind-effect corrections for 10V are similar to those for 6V.

Fig. 5 Conversion to SST from 6V or 10V

4. Algorithm Implementation
The AMSR2 SST algorithm is written with Fortran 77.

4-10

NDX-120015A

4.1.

Implementationcompile the source code by f77, and run “a.out”

4.2.

Input/output parameters
Input:

L1B
MGDSST / JMA for the initial SST
GANAL / JMA for air temperature

Output:SST

Quality flags for missing SSTs:
1 — land; land and coastal areas are removed
2 — sea ice; sea ice areas are removed for the input SST less than 5°C
3 — sun glitter; sun glitter areas of relative angle 25° or less with the sun are removed
4 — rain; rainy areas defined by “atmos_effect_6V” greater than 6.6 K are removed
5 — strong wind; areas of wind speeds greater than around 15 m/s are removed
6 — incidence angle; areas where the incidence angle deviates from 55° by 1° or more areremoved
7 — RFI; RFI areas from global stars and ground based radars are removed
8 — abnormal SST; SST less than –2°C

4.3.

Ancillary data
land flag (own data)climate salinity (own data)
scan bias for 6V and 6H (own data)atmospheric correction table (own data)monthly and SST dependency sp specification table (own data)brightness temperature correction table during orbit, among months, among years, if
necessary (own data)
average ASD table (own data)

4.4.

Processing flow
run several subroutines sequentially

4.5.

Example output
see the AMSR-E SST

4.6.

Limitations

4-11

NDX-120015A

range –2 to 35°C
no data for wind speed greater than around 15 m/s
no data for rain rate greater than about 0.5 mm/h
no data in areas within 100 km of the shoreline
no data in areas within 100 km of sea ice
no data inside sun glitter areas
no data within RFI areas
ground based radar
global stars

5. Validation Concept
We will compare the AMSR2 SST retrievals with buoy SST measurements observed for the
global ocean. Buoy data are collected through the Global Telecommunications System (GTS) andinclude data from both moored and drifting buoys. The number of matched observations is roughly
800 per day.
In order to compare the SSTs from AMSR2 and buoys, we will average the AMSR2 SST withineach 3 × 3 pixel area. We omit cases in which the difference between the maximum and minimum
AMSR2 SST is larger than 3°C within the 3 × 3 pixel area. We also omit cases where the absolutedifference between the AMSR2 and buoy SST is larger than 3°C.
We will check for cross-talk of the SST retrievals with four parameters: accumulated water
vapor content (WV), accumulated cloud liquid water content (CLW), sea surface wind speed, and
RWD.

References
Hollinger J. P. and Lo R. C., 1984: Determination of sea surface temperature with N-ROSS. Naval
Research Laboratory Memorandum Report 5375.

Milman A. S. and Wilheit T.T., 1985: Sea surface temperature from the Scanning Multichannel
Microwave Radiometer on Nimbus 7. J.G.R., 90 (C6): 11631-11641.

Klein L. A. and Swift C. T., 1977: An improved model for the dielectric constant of sea water at
microwave frequencies. IEEE Trans. Antennas and Propagation, 25: 104-111.

Shibata A., 1994: Determination of water vapor and liquid water content by an iterative method.
Meteorol. Atmos. Phys., 54: 173-181.

4-12

NDX-120015A

Shibata A., Imaoka K., Kachi M., and Murakami H., 1999: SST observation by TRMM Microwave
Imager aboard Tropical Rainfall Measuring Mission. Umi no Kenkyu, 8: 135-139 (in Japanese).

Shibata A., 2003: A change of microwave radiation from the ocean surface induced by air-sea
temperature difference, Radio Science, 38 (4): 8063-8072.

Shibata A., 2004: AMSR/AMSR-E SST algorithm developments: removal of ocean wind effect,
Italian J. Remote Sensing, 30/31: 131-142.

Shibata, A., 2006: Features of ocean microwave emission changed by wind at 6 GHz, J. Oceanogr.,62(3), 321-330.

Shibata, A., 2007: Effect of air-sea temperature difference on ocean microwave brightness
temperature estimated from AMSR, SeaWinds, and buoys, J. Oceanogr., 63(5), 863-872.

4-13

NDX-120015A

Chapter 5.
Description of GCOM-W1 AMSR2
Sea Surface Wind Speed Algorithm

Akira Shibata
Meteorological Research InstituteJapan Meteorological Agency

5-1

NDX-120015A

Table of Contents

1.

Introduction

3

2.

Algorithm Overview

3

3.

Theoretical Description

3

4.

Algorithm Implementation

85.

4.1 Implementation

8

4.2 Input/output parameters

8

4.3 Ancillary data

8

4.4 Processing flow

8

4.5 Example output

9

4.6 Limitations

9

Validation Concept9References

5-2

NDX-120015A

1. Introduction
There are two types of algorithm used for retrieving sea surface wind speed (SSW) from passive
microwave radiometers. One uses a higher frequency, such as 36 GHz, and the other uses lower
frequencies, such as 6 and 10 GHz. Although the first method has a better retrieval accuracy than thesecond method, it suffers significant contamination from raindrops. Consequently, the SSW retrievedby the first method is absent for rainy areas. An algorithm for the first method for the Advanced
Microwave Scanning Radiometer (AMSR) is described in Shibata (2012), and for the second methodin Shibata (2006). In developing the first type of algorithm, we used a combination of data fromAMSR and SeaWinds, both onboard the ADEOS-II. Here we present an algorithm using the higher
frequencies of AMSR2, onboard the GCOM-W1 satellite, which is the first satellite of Japan’s
Global Change Observation Mission.

2. Algorithm Overview
To develop the SSW algorithm, we need to know how the characteristics of brightness
temperature (Tb) vary with sea surface wind speed. Although it is more desirable to know this for afrequency of 36 GHz, it is difficult to determine the variation for such a high frequency, owing to theincreased contamination by atmospheric effects. Therefore, we determine it at 6 GHz. As shown in
Fig. 2 in the sea surface temperature (SST) algorithm description, the increment in Tb with verticalpolarization at 6 GHz (6V) because the wind speed is zero under weak wind conditions (≤ 6 m/s) andincreases gradually above this speed. The Tb with horizontal polarization at 6 GHz (6H) increases
even under weak wind conditions, and its increment is larger than that of 6V for all wind speeds. An
anisotropic feature depending on the relative wind direction (RWD), the angle made by the directionof AMSR viewing and the wind direction, is found for both polarizations. 6V reaches a maximum inthe upwind direction and a minimum in the downwind direction; 6H reaches a maximum in the
crosswind direction and a minimum in both the up and downwind directions. We assume that Tb
maintains similar characteristics at 36 GHz. Firstly, we will calculate a parameter, S36, from the 36
GHz data. Then, we will correct the anisotropic feature of S36, and finally, we will convert S36 to
SSW.

3. Theoretical Description
Calculation of S36

5-3

NDX-120015A

Tb at 36 GHz was computed at the satellite height using a microwave transfer model that has been
described in detail (Shibata, 1994; Shibata, 2004). The atmospheric part of this model deals onlywith absorption processes by water vapor, oxygen, and cloud liquid water, and does not deal with
scattering by raindrops. The reflection from a calm ocean surface is described by the Fresnel formula,for which the complex dielectric constant of ocean water is obtained from Klein and Swift (1977).
The model uses the daily aerological data of about sixteen Japanese observation sites in one year forintervening atmospheric parameters. Figure 1 reports the calculated Tb with vertical polarization at
36 GHz (36V; ordinate) versus Tb with horizontal polarization at 36 GHz (36H; abscissa), both inunits of Kelvin. The dots represent averaged values of 36H corresponding to 36V, with a 1 K interval.
In these calculations, the total amount of water vapor varies from 5 to 60 kg/m2, and the liquid watercontent takes two values of 0 and 0.4 kg/m2. The line OA is a regression line of these data. The line
OA represents the calm ocean surface condition, while the line WA represents the ocean surface
roughened by wind, for which the emissivity of the vertical polarization increases by 0.01 and that of
the horizontal polarization increases by 0.1. The two lines OA and WA intersect at the point A. For
other changes in surface emissivity, the regression lines also intersect at the point A. Tb at the point
A is about 300 K for both polarizations, which may represent the intervening atmospherictemperature in a lower level. Under an infinite nonopaque condition, Tb may approach theatmospheric temperature in the lower level, regardless of the ocean surface condition.

Fig. 1 Definition of S36

A parameter representing ocean wind is defined by a horizontal length between OA and WA, as

5-4

NDX-120015A

shown by S36 in Figure 1, and written as in eq. (1).
S36 = (36H – a × (36V – b) – c) / f + t(1)

b = 208
f = 1 – 0.01 × (36V – 200)

(2)

The coefficient b is a constant, whereas a and c depend on SST. In eq. (1), f is an atmosphericcorrection and represents the convergence of the two lines OA and WA. The difference in Tb for 36V
between O and A is about 100 K, and f diminishes when 36V approaches 300 K in eq. (2). In eq. (1),t is an offset that will be explained later. The position of point O shifts in accordance with changing
SST, and the values of a and c depend on SST. The values of a and c can be determined bycomparing S36 with the SeaWinds wind speed. Ideally, these values should be determined against a
SeaWinds wind speed of 0 m/s, in which case t = 0. However, such a case is very rare, and there areinsufficient collocated data points for comparison. Therefore, collocated data are obtained at a
SeaWinds wind speed of 5 m/s (i.e., in a range from 4.5 to 5.5 m/s), and t is set empirically as 4.5.
Table 1 lists values of a and c for seven SSTs from 0 to 30ºC with 5ºC intervals. The techniqueemployed in eq. (1) can be applied also to the other frequencies of 18, 23, and 89 GHz, although
their SST dependencies are not listed. Hereafter, the corresponding parameters for these frequencies
are called S18, S23, and S89, respectively.

Table 1 - Values of a and c against seven SSTs
SST [ºC]

a

c

0

2.23

132.0

5

2.20

132.2

10

2.14

131.5

15

2.07

130.7

20

2.06

128.8

25

2.03

127.4

30

2.06

124.2

Correction of the anisotropic feature of S36
S36 and the corresponding parameters for the other frequencies also have anisotropic features
depending on RWD, whose dependencies can be determined by comparing the parameters with the
SeaWinds wind speed and direction. In Figure 2, S89, S36, S23, and S18 are shown against RWD
and for ten cases of SeaWinds wind speed from 1 to 19 m/s with 2 m/s intervals. It is clear that thesensitivity to wind speed becomes greater with increasing frequency, i.e., highest for S89, medium
for S36, and lowest for S23 and S18. However, the behavior of S89 seems to be noisier than for the

5-5

NDX-120015A

other frequencies, which may be attributable to atmospheric contamination at higher frequencies.
Therefore, S36 is the best parameter for retrieving wind speed. S36 is almost constant with RWD for
weak wind speeds of less than 5–7 m/s. For wind speeds over this value, S36 has an anisotropic
feature, reaching a minimum in the upwind direction and a maximum in the downwind direction. As
mentioned above, 6V reaches a maximum in the upwind direction and a minimum in the downwind
direction. Assuming that this feature of 6V holds similarly for 36V, S36 reaches a minimum in the
upwind direction, because a × 36V is dominant in eq.(1) (a is almost equal to 2). As for the
maximum value, S36 reaches a maximum in the downwind direction, and S18 reaches its maximum
for directions between cross and downwind directions (S23 has characteristics between those of S36
and S18).

Fig. 2 S89, S36, S23, and S18 against RWD and for ten cases of SeaWinds wind speed.

The anisotropic feature of S36 should be corrected before converting S36 to wind speed, because
otherwise the converted wind speed will inherit the anisotropic feature. Figure 3 depicts the relation
between S36 and 6H for three cases of RWDs (up, cross, and downwind) and ten cases of windspeeds (from 1 to 19 m/s with 2 m/s intervals). 6H is defined by eq. (2) in the SST algorithm
description. As mentioned previously, S36 reaches a maximum in the downwind direction and aminimum in the upwind direction. Also, 6H reaches a maximum in the crosswind direction and aminimum in the up and downwind directions. The anisotropic correction means that S36 for both theup and downwind directions is adjusted to S36 for the crosswind direction, following the slopes in

5-6

NDX-120015A

Figure 3. A fitting curve corresponding to crosswind is shown in Figure 3. Different slopes are used
for the adjustments from up to crosswind and from down to crosswind directions, respectively.
Figure 4 depicts S36 before and after the anisotropic correction. For almost the entire range of
SeaWinds wind speed, S36 remains constant with RWD. However, at larger wind speeds, such as 19
m/s, S36 varies slightly with RWD, with variations on the order of 1 K.

Fig. 3 Method of correcting the anisotropic feature of S36

Fig. 4 Comparison of S36 before and after anisotropic correction

Conversion of S36 to wind speed
Finally, S36 is converted to wind speed, using the curve shown in Figure 5, which was obtainedby comparing S36 with the SeaWinds wind speed for about 3 × 108 collocated data points over a

5-7

NDX-120015A

period of seven months. Conversion of negative S36 values is considered, as S36 sometimes
becomes negative. Reasons for S36 to become negative are (a) Tb noise (fluctuation) of the AMSR
and (b) incorrect SST derived from the weekly SST. The remaining error in the corrected S36 (Fig.
4) gives rise to a wind speed error. From the conversion curve, an error of 1 K in S36 is estimated as
an error of 1 m/s in wind speed in the 15–20 m/s range.

Fig. 5 Conversion curve from S36 to wind speed.

4. Algorithm Implementation
The AMSR2 SSW algorithm is written in Fortran 77.

4.1.

Implementationcompile the source code with f77, and run “a.out”

4.2.

Input/output parameters
Input:

L1B
MGDSST / JMA for the initial SST
GANAL / JMA for air temperature

Output: SSW
Quality flags for missing SSW:
1 — land; land and coastal areas are removed

5-8

NDX-120015A

2 — sea ice; sea ice areas are removed for input SST less than 5°C
3 — sun glitter; sun glitter areas of relative angle 25° or less with the sun are removed
4 — rain; rainy areas defined by “atmos_effect_6V” greater than 6.6 K are removed
5 — abnormal wind; impossible to determine the relative wind direction
6 — incidence angle; areas where the incidence angle deviates by 1° or more from 55° are
removed
7 — abnormal wind; out of range from 0 to 30 m/s
8 — not used

4.3.

Ancillary data
land flag (own data)
S36 specification table (own data)
conversion table (own data)
Tb correction table during orbit, among months, among years, if necessary (own data)

4.4.

Processing flow
run several subroutines sequentially

4.5.

Example output
see the AMSR-E SSW

4.6.

Limitations
range of 0 to 30 m/s
no data for rain rate greater than about 0.5 mm/h
no data for areas within 20 km of shorelines
no data for areas within 20 km of sea ice areas
no data inside sun glitter areas

5. Validation Concept
We will validate the AMSR2 SSW by comparison with SSW values from ocean buoys. Ocean buoy
data are collected from moored buoys of the National Data Buoy Center (NDBC), the TropicalAtmosphere Ocean (TAO) Project, and the Prediction and Research Moored Array in the Atlantic(PIRATA) in open ocean 100 km off shorelines. The height at which the wind speed is measured
differs among buoys, so it is adjusted to 10 m for all buoys by assuming a neutral atmosphericcondition. In combining the AMSR2 and buoy data, the time difference between the AMSR2 andbuoy measurements is limited to within one hour. Nine pixels of AMSR2 data are averaged aroundthe buoy location.

References

5-9

NDX-120015A

Klein L. A. and Swift C. T., 1977: An improved model for the dielectric constant of sea water at
microwave frequencies. IEEE Trans. Antennas and Propagation, 25: 104-111.
Shibata A., 1994: Determination of water vapor and liquid water content by an iterative method.
Meteorol. Atmos. Phys., 54: 173-181.
Shibata A., 2004: AMSR/AMSR-E SST algorithm developments: removal of ocean wind effect.
Italian J. Remote Sensing, 30/31: 131-142.
Shibata, A., 2006: A wind speed retrieval algorithm by combining 6 and 10 GHz data from
Advanced Microwave Scanning Radiometer: wind speed inside hurricanes. J. Oceanogr., 52(3),pp.351-359.
Shibata, A., 2012: Ocean wind speed retrieval algorithm using the frequency 36 GHz
vertical/horizontal and 6 GHz horizontal data of the Advanced Microwave Scanning Radiometer(AMSR). Italian J. Remote Sensing, 45: 133-140.

5-10

NDX-120015A

Chapter 6.
Description of GCOM-W1 AMSR2
Sea Ice Concentration Algorithm

Josefino C. Comiso
Goddard Space Flight Center
National Aeronautics and Space Administration

Kohei Cho
Tokai University
6-1

NDX-120015A

Table of Content

1.

Introduction ................................................................................................................... 3

2.

Algorithm Overview ....................................................................................................... 4

3.

Theoretical Description .................................................................................................. 5

4.

Algorithm Implementation ........................................................................................... 20

5.

4.1.

Implementation ....................................................................................................... 20

4.2.

Input/Output parameters .......................................................................................... 21

4.3.

Ancillary data ......................................................................................................... 21

4.4.

Processing flow ...................................................................................................... 21

4.5.

Example output....................................................................................................... 24

4.6.

Limitations ............................................................................................................. 24

Validation Concept ....................................................................................................... 24

References ........................................................................................................................... 25

6-2

NDX-120015A

1. Introduction
Sea ice is an important component of the Earth's climate system.

At any one time, sea ice covers

about 6% of the world’s oceans which are the primary reservoir of global heat and energy.
influences the oceans and the atmosphere in many ways.

Sea ice

In particular, the growth of sea ice in thepolar regions enables the formation of high density bottom water that becomes part of the global
thermohaline ocean circulation (Broecker, 1991).

The sea ice cover also changes the spatial

distribution of salt in the oceans in that where it forms, the salinity of the water increases because ofbrine rejection and where it melts, the salinity of the water decreases because of the introduction of
low salinity melt water (Wadhams, 2000).

The low salinity melt water at the same time enables theformation of a stable layer at the surface exposed to abundant sunlight, and with ample nutrients andiron can be the ideal platform for photosynthesis in marine plants (Walker and Comiso; 2009).

Thisphenomenon is in part the reason why the high latitude regions are among the most biologically
productive regions in the world.

Sea ice is also a good insulator, keeping the ocean heat from being

released to the atmosphere in autumn and winter and also keeping the ocean from absorbing too
much heat in the spring and summer (Weeks and Ackley, 1986).

It is this phenomenon, called

ice-albedo feedback, that causes the amplification of global climate change amplified in the polar
region and that makes sea ice the barometer for the global climate.
Historical records of the global sea ice cover, as derived from satellite passive microwave
data, have provided some of the most important climate change signals that may be associated withanthropogenic greenhouse gases.

One of the most revealing results from satellite microwave datais the large rate of decline of the Arctic perennial ice extent which is currently at -13.5 % per decadeand its thick multiyear ice component at -15% per decade (Comiso, 2012).

In September 2007, theperennial ice cover was not only the record low but also 40% lower than the average value since1978, and 27% lower than the previous low value in 2005 (Comiso et al., 2008). A slight recovery
followed but in September 2011, the perennial ice cover was about as low as in 2007. Thisphenomenon has been recognized as one of the most intriguing climate change signals coming fromthe Arctic region.

The anomalously high sea surface temperature during the period and the strong

possibility of having an ice free ocean in the Arctic in the summer in the near future has beenreported (Shibata et al., 2009; Lindsay and Zhang, 2005; Wang and Overland, 2009).

The potential

impact of such an event on the ecology, environment, navigation and mineral exploration has been
the subject of many studies (e.g., Bhatt et al., 2010; Smith and Comiso, 2009). On the other hand,similarly intriguing is that the sea ice cover in the Antarctic region has been changing as well but inthe opposite direction (Zwally et al., 2002; Comiso et al., 2009).

Such positive trend, which is now

about 2% per decade, has been postulated to be in part caused by the deepening of the lows in West
Antarctica associated with the Ozone hole that in turn have caused stronger winds off the Ross ice
shelf and hence more ice production in the Ross Sea region (Turner et al., 2009).

6-3

Enhanced ice

NDX-120015A

production in the Ross Sea has been confirmed (Comiso et al., 2009) but other studies indicate that
other processes not linked to stratospheric ozone depletion must be invoked to explain the observed
increase in the Antarctic extent since 1979 (Sigmond and Fyfe (2010).

Because of many important

issues associated with the aforementioned changes, accurate characterization of the global sea ice
cover more than ever needed. Such product should also be consistent with historical data to enable
meaningful trend and variability studies.

2. Algorithm Overview
The Bootstrap Algorithm has been used for generating the standard sea ice concentration
products from AMSR-E and now for AMSR-2 data.

AMSR-2 data will provide continuity tocurrent AMSR-E data that have become the baseline and primary source of sea ice cover. The
algorithm description will cover the following: (a) the physical and theoretical basis for the estimateof ice concentration; (b) details about the determination of the reference brightness temperatures for
100% and 0% sea ice concentrations; and (c) the techniques used for land and open ocean masks;
and show (d) sources of errors and the consistency of the data with historical AMSR-E and SSM/I
data.

Potential refinements in the algorithms associated with better resolution, improved calibration,less radiometric noise and more frequency channels will also be presented.
Monitoring the sea ice cover with satellite passive microwave data has had a long history,
starting with the Nimbus-5 Electrically Scanning Microwave Radiometer (ESMR) launched in
December 1972, continuing with the Nimbus-7 Scanning Multichannel Microwave Radiometer(SMMR), launched in October 1978, and then with the series of DMSP Special Scanning Microwave
Imagers (SSM/I) the first of which was launched in July 1987.

The data from these sensors haveprovided the continuous and consistently derived sea ice concentration data from 1978 up to thepresent (Comiso and Nishio, 2008).

The data, however, are not without any shortcomings among

which are the relatively coarse resolution, significant radiometer noise and antenna side lobes.

Thelaunch of the Japanese Advanced Microwave Scanning Radiometer (AMSR-E) on board theEOS/Aqua satellite in May 2002 and a similar version (AMSR) on board ADEOS-2 satellite in
December 2002 provided a new era in the remote sensing of sea ice.

The new system has a muchbetter resolution, wider spectral range and larger swath width than previously available data.

For

example, the resolution of the 37 GHz channel for SSM/I is 30 x 38 km while that of AMSR-E atapproximately the same frequency is 10.3 x 13.7 km.

The higher resolution has enabled moreaccurate identification of the ice edge and improved characterization of the ice cover in the innerzones (Comiso and Nishio, 2008).

It also enabled assessments of snow cover and surface ice

temperature as well as improved estimates of drift velocities of the ice floes (Markus and Cavalieri,
1998; Comiso et al., 2003; Yaguchi and Cho, 2009).

Further improvements in accuracy are

expected with the launch of the GCOM-W/AMSR2 sensor because of significantly higher resolution

6-4

NDX-120015Aof the latter (i.e., 7x12 km for AMSR2 compared to 10x14 for AMSR-E and 30x38 km for SSM/I at
37 GHz).
The Bootstrap Algorithm has had a long history and has been used to generate sea ice datafrom 1978 to the present.

It makes use of three frequency channels (i.e., 18.7 at verticalpolarization, 36.5 GHz at both vertical and horizontal polarization and 22 GHz at vertical resolution)that are available continuously from SMMR, SSM/I, AMSR-E.

The launched of AMSR2 aboard

the GCOM-W satellite in May 2012 is a most welcomed development, especially with the demise of
AMSR-E in 2010. The algorithm is thus ideally suited for generating long term and consistent data
set needed for detailed studies of the variability and trends in the sea ice cover.

Among theobjectives of the Global Change Observation Mission (GCOM) seeks to understand global climatechange through the use of long-term satellite observing systems.

It is important that the sea ice data

set that had been enhanced by AMSR-E will continue with AMSR2, especially during the current
climatic period when the sea ice cover is undergoing not just large seasonal and interannual
variability but also rapid declines in the Arctic and some increases in the Antarctic.

3. Theoretical Description
The conversion of digital satellite data to sea ice parameters requires a knowledge of somebasic geophysical parameters, among which are surface temperature and the microwave emissivityof the surface at the different frequency and polarization channels of the sensor.

In addition,

atmospheric effects that vary from one location to another and from one day to the next have to betaken into account. The brightness temperature, (T B ), recorded by satellite passive microwave
sensors at a given wavelength can be estimated using the basic radiative transfer equation given by
τ

τ

0

0

TB = εTSe −τ + ∫ T(z)ς (z)e −τ +τ'(z)dτ '(z) + (1 − ε )κe −τ ∫ T(z)ς (z)e −τ '(z)dτ '(z)(1)

where ε is the emissivity of the surface, T S is the physical temperature of the surface, τ’(z) and τ arethe atmospheric opacities from the surface to a height z and from the surface to the satellite height,
respectively, κ is an estimate of the diffusiveness of the surface reflection, and T(z) and ζ (z) are thetemperature and the emittance at z. In equation (1), the first term represents radiation directly fromthe earth’s surface, which is the dominant contribution for measurements at microwave frequencies.
The second term represents satellite observed upwelling radiation that emanates directly from theatmosphere, while the third term represents downwelling radiation from the atmosphere that has beenreflected by the surface before it reaches the satellite detector. A fourth term that takes into accountthe reflected contribution of radiation from free space (i.e., the cosmological 2.7 K contribution fromthe Big Bang), which is an additive contribution, is usually negligible and not included in equation (1).

6-5

NDX-120015A

This radiative transfer equation provides the means to transform the top of the atmosphere brightness
temperature as measured directly by satellite sensors to surface brightness temperature needed toobtain the geophysical surface parameter. The brightness temperature of the surface of interest isequal to εT S which is part of the first term in equation (1). The atmospheric parameters and alsosurface emissivity as required by the equation changes spatially and temporally and are unfortunately
not known on a real time basis and hence the need to develop an alternative procedure. The Bootstrap
Algorithm is designed to take the atmospheric contribution into account in an indirect manner througha technique based on analysis of real time data. With the Bootstrap Algorithm, a multichannel
scheme was developed in which the patterns in the distribution of data points and clustering of thesepoints in multi-dimensional space are used to indirectly infer emissivity and at the same time, take into
account spatial variations in surface temperature and atmospheric effects.

2.3.1 Tie-Point and Threshold Optimization
The basic assumption of the Bootstrap Algorithm for ice concentration is that within theice pack, the surface is covered by either sea ice or ice free (liquid) water.

The brightness

temperature recorded by the satellite passive microwave sensor within each field-of-view is thus acontribution of radiation from ice covered surfaces and from ice-free surfaces. The brightness
temperature, T B , as measured by AMSR-2 is thus assumed to be a linear combination of the two
surfaces as expressed by the following mixing formulation:

TB = TICI + TOCO

(2)

where T I is the brightness temperature of 100% ice covered areas and T O is the corresponding valueof 100% open water, while C I and C O are the percentage concentrations for sea ice and open (icefree) water respectively. Using the relationship, C O + C I= 1, and assuming that the observational

area is covered only by either ice or water, equation (2) can be solved explicitly for C I and the result
is as follows

CI =TB − To
TI − TO(3)

Equation (3) is the fundamental formula that has been used by all sea ice concentration algorithms.
It is a relatively simple formula but not so easy to implement.

The big challenge has been how to

obtain good estimates of the surface values of T B , T O and T I which are all functions of emissivity(ε), surface temperature (T S ), and atmospheric opacity.

In principle, the equation is independent ofsensor resolution and should provide similar results regardless of resolution.

6-6

But there are

NDX-120015A

exceptions especially near boundaries of ice and water and in areas where the contrast in theemissivity of sea ice and water are not strong enough to enable a good discrimination of the two
surfaces.

In such areas, there could be a smearing of data as is normally the case at the marginal ice

zone in part due to antenna side lobes.

In general, the higher the resolution, the better the chancesin getting the correct location of the ice edge and in getting accurate characterizations of themarginal ice zones.
An ideal algorithm for retrieving the ice concentration would be the one that calculates
accurately the surface values of the parameters T B , T O and T I in equation (3) at each field-of-view(or data element) of satellite observation.

This means that all three parameters have to becorrected for atmospheric effects at each measurements.

The usual technique is to utilize aradiative transfer model which makes use as input, atmospheric profile parameters to take into
account the effect of spatially and temporally changing opacity τ of the atmosphere as indicated inequation (1). Such atmospheric contributions can be as large as 20% of the satellite observed
signal and are therefore important to take into account.
Earlier algorithms were based on radiative transfer models that used atmospheric profiles
taken from averages derived from radiosonde data at some Arctic regions (Svendsen et al., 1983;
Swift et al., 1985).

These algorithms were reasonably successful in regions where and time periods

when such radiosonde data are acquired and used for atmospheric corrections.

Generally, however,there are problems when they were used as a global algorithm because of large spatial and temporalvariabilities of atmospheric and surface conditions.

The algorithms that are currently more

frequently utilized are those that take advantage of the multichannel capability of the satellite data in
obtaining the required input parameters in equation (3) as described in Comiso et al. (2003). Twoof these algorithms are the Bootstrap Algorithm and the Nimbus-7 Team Algorithm (now called theNASA team algorithm), both developed at the NASA Goddard Space Flight Center (Cavalieri et al.,
1984; Comiso, 1986).

The original Team Algorithm used three tie points: two for the two general

ice types (i.e., first year and multiyear ice) and one for open water, and employs polarization andgradient ratios to minimize the effect of varying surface temperatures. The different techniques foraccounting for spatial changes in ice temperature and emissivity and the use of different sets ofchannels yielded different results between the Team Algorithm and the Bootstrap Algorithm (Comiso
et al., 1997; Comiso and Steffen, 2001) in large areas of the polar regions.

The different results

were attributed to different sensitivity of the polarization ratios to surface effects, including ice andsnow layering, as discussed in Matzler (1984).

The original NASA Team Algorithm (now calledNT) has been substantially revised to make use of the 89 GHz in combination with other channels
and is currently called NT2 (Markus and Cavalieri, 2000).

With NT2, a radiative transfer model(e.g., Kumerow, 1993) has to be used to correct for the high sensitivity of the 89 GHz data toatmospheric effects.

It also requires some accounting on the large variations in the surface

6-7

NDX-120015A

emissivity of sea ice at this frequency.

The results from NT2 have been shown to be an

improvement to that of NT and are in much better agreement with those of the Bootstrap Algorithm
but there are still significant differences (Parkinson and Comiso, 2008; Comiso and Parkinson, 2008).
Such differences are believed to be associated in part with the use of the 89 GHz data (in NT2)
which at times could cause errors due to inadequate accounting of atmospheric and surface effects.
The differences are also in part due to differences in resolution since the two algorithms uses
different sets of frequency channels.
The Bootstrap Algorithm takes advantage of the unique distribution of brightness
temperature (or emissivity) data points in two or three dimensional space to assess and evaluate
patterns that may be associated with the general distribution of data points in sea ice covered regions.
An schematic illustration of such a distribution is shown in Figure 1 using data from two radiometerchannels. Data points in consolidated ice regions where the ice concentrations are close to 100% areusually form a linear cluster along a line near AD.
the line AW.

Those in the open water areas would be along

The scheme for finding the two tie points for ice concentration is to utilize the line

AD in Figure 1 as the reference for 100% ice and the data point O for 100% liquid (or open) water.
Data points located between the line AD and O are interpreted as having ice concentrations between
0 and 100% in this scheme.

Thus, given a data element at a hypothetical data point B in the plotmeasured by the satellite sensor, the ice concentration can be derived by first extending the line
along OB until it intersects the line AD. The brightness temperature T I at the intersection point, I,
represents 100% ice for this particular ice type/surface while the brightness temperature T O at O
represents 0% ice. Using this tie-point values in equation (3), ice concentration can be calculated.
This can be done for either channels (along the horizontal or the vertical), but when T O and T I

areclose or equal to each other, as can happen at 36 GHz (V) and higher frequencies, there is asingularity problem (i.e., if T O = T I ). Thus, the ratio of OB to OI is usually used, as described inComiso (1995), which provides the same value as the ratio of the numerator and denominator inequation (3). It should be noted that ice concentrations below a cut-off of about 10% are derived
because the emissivities of ice and water are mainly impossible to discriminate at such values as
discussed below.

6-8

NDX-120015A

Figure 1.

Schematics for the choice of tie-points in the Bootstrap Algorithm.

Data points along

AD are considered as having 100% ice concentration while those along OW have 0% ice
concentration.

The line AD is determined by the algorithm on a daily basis and inferred from a regression
analysis of data points above the red dash line in Figure 1.

The red dash line is determined suchthat most of the data points above it have near 100% sea ice concentration.

A small positive valueof a few Kelvins is added to the offset of the regression line to account for the known presence ofopen water (of about 2 % in the Arctic and 4% in the Antarctic) for much of the consolidated ice
regions.

The estimate of the parameters of the line AD is critical to the accuracy of the ice

concentration estimate.

Once AMSR2 data are available, the values of the offsets will be optimized

such that the retrieved ice concentration values are as consistent as possible with those fromvalidation data sets (e.g., high resolution satellite data or detailed ship or aircraft observations) and at
the same time allows for variations in surface emissivity and atmospheric attenuation.
The criteria used for choosing the channels for the algorithm are: (a) they provide thegeophysical parameter at the optimal resolution and accuracy; and (b) ancillary data are used only if
absolutely necessary to improve accuracy.

Thus, the set of 36.5 GHz channels (called HV36) isused because of reasonably good resolution and at the same time takes into account spatial variations
in surface temperature and emissivity for accurate retrieval. With this set, the data in the
consolidated ice regions form a linear cluster with a well defined slope for AD (when one

6-9

NDX-120015A

polarization is plotted versus the other) the value of which is approximately equal to one. The
values form a linear distribution because the wavelengths are the same and emissivity of ice isapproximately the same for the two polarizations.

Thus, they are affected by intermediary factors,

like snow and the atmosphere, in much the same way.

What makes it even more valuable is thatthe set is basically insensitive to spatial variations in temperature since the slope is approximately
one and any change in th brightness temperature due to temperature as observed in one channel isapproximately equal to that in the other channel.

Thus the net effect of a changing temperature is tocause T I to slide along the line AD and practically no impact on the accuracy in the retrieval of ice
concentration.
The algorithm complements the use of the set of 36.5 GHz channels with another set using
18 (or 19) GHz at vertical polarization which we call th V1836 set.

This other set provides themeans to remove some ambiguities in the HV36 set since in this set of channels, the open ocean datacluster (i.e., the line OW) is either along the line OA in Figure 1 or to the left of this line.
Furthermore, the horizontal channel is more sensitive to layering and other surface effects than thevertical channel (Matzler et al., 1984) and in some ice covered areas, the data points fall below or to
the right of the line AD.

Both problems are resolved through the additional use of a set of channels

that utilizes the 18 (or 19) GHz in combination with the 36.5 (or 37) GHz channel at verticalpolarization (thereafter called the V1836 set) which provides an even easier discrimination and agood contrast between the emissivity of ice and water.

The use of the V1836 set shows somesensitivity to variations in surface ice temperature but the associated error is relatively minor because
sea ice is usually covered by snow which is a good insulator and the spatial variability of observed
surface ice temperature is not much.

The standard deviation of observed snow/ice interface

temperatures has been about 2.5K.

2.3.2 Statistical Error Analysis
Figure 2a provide the means to evaluate how accurately the sea ice concentration can be derived
from the HV36 set of channels.

The data points are from winter data in the Northern Hemisphere andthe distribution for consolidated ice (in blue) is indeed quite compact and the slope is close to 1. The
accuracy in the estimate of C I depends mainly on the accuracy in the estimate of T I

since it is known

that T I varies a lot more than T O in the pack ice regions. The accuracy in the estimate of T I
dependent on how well AD represents 100% ice.

is in turn

The more well defined the linear cluster AD is, themore accurate the retrieval is going to be. To quantify how well defined the ice cluster AD is, the scatterplot is rotated such that the AD cluster is along the vertical as shown in orange in Figure 2a. The widthof the cluster can then be quantified with a frequency histogram of the sum of data points along thevertical within each horizontal bin and the results are shown in Figure 2b. In the plot, the peak in the left
represents the distribution of consolidated ice data points along an arbitrary horizontal axis that can be

6-10

NDX-120015A

converted to ice concentration since the relative location of TO is also shown in the plot (approximately
the highest data point to the right). After the conversion, the standard deviation of the ice peak is about
1.5 units which is estimated to be equivalent to less than 3% ice concentration. The actual uncertainty inthe ice concentration associated with the variability of

TI

would be even smaller since theconsolidated ice data points include a fraction of open water within a relatively large footprint (which inthis case is 25 by 25 km).

Figure 2. (a) Scatter plot of 36 GHz(V) versus 36 GHz(H) brightness temperatures with the ice cluster
along AD shown as blue data points. Also shown is the same set of data points with data points rotated(in orange) such that the line AD cluster is along the vertical.
points (orange) in (a).

(b) frequency distribution of rotated data

The standard deviation of the distribution for the ice curve (red line) is 1.5 oK

which is approximately 2% in ice concentration.

2.3.3 Systematic and Residual Errors

6-11

NDX-120015A

2.3.3.1

Ocean and Land Masks
The aforementioned retrieval technique for ice concentration is appropriate only in ice

covered oceans.

Applying the same algorithm over land and open oceans would yield not zero ice

concentration values and hence the need for a land and liquid ocean masks.

Such land and liquidmasks are usually not perfect and would affect the accuracy of the retrieval of areas covered by sea
ice.

The land mask, that has been previously used were based on published boundaries of land

areas. The published boundaries have been shown to be inaccurate in some areas when comparedwith high resolution satellite data.

This is in part because in the polar regions, the land/ocean

boundaries are not permanent because of iceberg calving, glacial surging or retreat and ice melt.
The land mask to be used for AMSR2 will be an enhanced version of the one currently used for
AMSR-E and will be based on MODIS or higher resolution satellite data.
An open ocean mask has been developed by identifying the location of data points that
correspond to such open ocean areas in the scatter plots of one channel versus another.

Most of thedata points in the open ocean can be isolated and identified readily using this procedure but there are
significant number that are problematic. In particular, areas in the open ocean that are under the
influence of inclement weather conditions can have signatures similar to those of ice covered ocean.
We make use of a combination of 19, 22, and 36 GHz brightness temperatures (T B ) at vertical (V)and horizontal (H) polarizations to discriminate open ocean data from ice data during such unusual
conditions as discussed in Comiso (1995).

Figures 3a and 3b show scatter plots of T B (19,V) versus

the difference T B (22,V) – T B (19,V) using SSM/I and AMSR data, respectively.

The blue data

points in the scatter plot along a line OW represent data from the open ocean at all weather
conditions while the black data points are those from ice covered ocean. In the open ocean, strong
winds and inclement weather cause large surface waves and foam that alters the emissivity of thesurface.

The microwave signature of open water is thus variable and moves from low values at O,

during relatively calm conditions, to higher values and towards W, in the scatter plot, depending on
the strength of the disruption of the surface water.
masked since they represent open water.

In the algorithm, these data points in blue are

The brightness temperature of open water within the ice

pack correspond to the lowest values of the cluster along OW since inside the pack the water surface
is relatively still. The tie point used to represent T O in equation (1) is thus a low value close to thelabel O (along OW).

The primary mask for open water is the set of channels shown in Figures 3aand 3b. This is usually complemented by the use of the V1936 set to remove residuals.

The sloping

red line represents approximately 10% ice concentration and is used as a threshold.

Below thisthreshold, ice and open ocean are not possible to discriminate as indicated by the data.

The vertical

red line is used to mask out data to the right of the line which represent primarily ice free water
under abnormal or extreme surface conditions, as described in Comiso et al. [2003].

We also use

climatological surface temperature data or SST from AMSR2 (Shibata et al., 2009) for unusually

6-12

NDX-120015A

difficult cases away from the ice pack assuming that the surface temperature of ice covered surfaces
cannot have surface temperatures greater than 5 oC.

Because of significant variability of theemissivity of sea ice near the ice edge, the error in the retrievals near 10% ice concentration isrelatively large compared to high ice concentration areas.

In the estimates for ice extent, which is

the integral sum of the area of data elements with ice concentration greater than a certain threshold.
We use 15% ice concentration as the threshold as was done in Zwally et al. [1983] (instead of 10%)for good consistency in the estimate of extents.

Among the refinements planed for the AMSR-2

algorithm is to set the threshold for ice extent at around 10% or as low as possible.

Figure 3.

Scatter plots of T B (V19, V) versus T B (22, V) - T B (19, V) for (a) SSM/I and (b)AMSR-E data. Data points in blue correspond to data interpreted as either ice free or <10% ice
concentration.

The sea ice algorithm is designed to work primarily in sea ice covered areas but since it is a
global algorithm, the algorithm is applied for all surfaces.

It is important to note that in theprocess, non-zero ice concentration are retrieved over land/ocean boundary areas far away from icecovered areas and where sea ice is not expected to be present.

At these boundaries, the satellite

observations of what is regarded as ocean area may be contaminated by emissions from land areas.
This has been a problem with SSM/I data and to a lesser extent with AMSR-E data because of thehigher resolution of the latter.

It is expected that with even higher resolution, AMSR2 will provide

even better accuracy in the region.

Data elements slightly farther away from the ice edge are also

6-13

NDX-120015A

affected because of antenna side-lobes which cause the brightness temperatures of data near the iceedge to be different when the satellite approaches the boundary from land compared to the otherdirection.

It is even more of a problem when there is no ice cover at the boundary because of themuch larger contrast between the emissivity of open ocean and land than between sea ice and land.
Thus, along the coast of places like Spain, there would be one or more pixels from the coast linesthat are retrieved as ice covered areas. A correction algorithm originated by Cho et al. (1996) was
adapted to mask out these problem areas.

The procedure is quite effective but not 100% effective,

especially where the spill over of land effects over the ocean areas are significant. The higher
resolution of AMSR2 data will further reduce the errors in these regions.
Only about 15% of the winter ice cover in the Southern Hemisphere survives the summer
to become perennial ice but even the latter apparently have signatures similar to those of seasonal ice(Zwally et al., 1983; Gloersen et al., 1992).

This is in part because the Antarctic perennial ice

consists mainly of second year ice which appears to have signatures closer to that of first year ice
than those of the multiyear ice in the central Arctic (Comiso, 2006).

The earlier version of the

Bootstrap algorithm took advantage of this and for simplicity and convenience, utilized only one setof channels (i.e., the V1836 set) in the retrieval of ice concentration.

Comparative studies withhigh resolution satellite visible and infrared data indicated that the retrieved values using the originalalgorithm has a negative bias of a few percentage ice concentration. Such bias was removed when
the algorithm used for the Arctic is also used in the Antarctic.

For optimum accuracy, however, thetie points used for Antarctic sea ice is determined independently from those used for Arctic sea ice.
Data from the AMSR-E on board the EOS/Aqua satellite which was launched in May 2002
provided many opportunities to validate the performance of the Bootstrap Algorithm that will be
used for processing AMSR2 data.

The most important validation tools have been the MODIS 250

m visible and near infrared data, and the AMSR-E 89 GHz data at vertical and horizontal
polarizations.

Figure 4 shows a side by side comparison of the sea ice cover in the Weddell Sea on

8 October 2004 as recorded by the 0.65 µm MODIS and ice concentrations from the AMSR-E 89
GHz channels at 6.25 km resolution and from the standard AMSR-E algorithm gridded at 12.5 km
resolution.

The 89 GHz ice concentration map was derived using the same Bootstrap ice

concentration algorithm technique but using the 89 GHz data at both horizontal and verticallypolarizations .

It is apparent that the ice cover as depicted by the 89 GHz image show spatiallydetailed features similar to those shown in the 250 m resolution MODIS data.

The three areas ofspring polynyas (lower middle) and the large leads in between are depicted in a similar fashion by
Figures 4a and 4b.

The marginal ice zone (top area) and the ice edge are also similarly represented

in the first two images except in areas where clouds are apparent in the MODIS data.

The standard

ice concentration map as depicted by Figure 4c shows good consistency with both MODIS and 89
GHz AMSR-E data but there are significant differences.

6-14

Because of the coarser resolution the

NDX-120015A

standard ice map does not show the lead patterns as distinctly as the other two and the spatial detailsof the polynyas and the marginal ice zones are not as well represented.

Also, some of the thinner

ice types, especially near the marginal ice zone, are represented as having intermediate ice
concentration values.

This is in part because the emissivity of thin ice is usually between that of

the thicker ice types and open water.

Furthermore, large waves near the marginal ice zone usuallycause the ice floes to be flooded and wet at the surface thereby causing alteration of the surface
emissivity.

However, the ice concentration map from the standard algorithms shows goodconsistency overall with those from MODIS and 89 GHz data. Among the refinements for the
AMSR2 Bootstrap Algorithm would be to find the optimum tie points that yields ice concentrations
that best match the MODIS and 89 GHz data.

Adjustments will also be made to make AMSR2 dataconsistent with that of AMSR-E sea ice data.

Figure 4.

Sea ice cover in the Weddell Sea on 8 October 2004 as represented by (a) MODIS 250 m

resolution data at 0.65 um; (b) ice concentration map from AMSR-E using the 89 GHz data athorizontal and vertical polarizations; and (c) sea ice concentration color-coded map using the
AMSR-E standard Boostrap Algorithm.

2.3.3.2

New Ice Areas and Pancakes

Another serious source of ambiguity is the presence of new ice, the emissivity of which varies withthickness and surface condition. Even in the simple case of an area covered mainly by consolidated
ice but with recently formed leads in mid winter, there is already a problem. If the fraction of lead
area is, for example, 10%, the estimated ice concentration should be 90%. This may indeed be thecase if the satellite passes by the area at the time when the leads were forming and the emissivity of thelead is the same as that of open water. The emissivity, however, changes with time as the same lead islikely covered by grease ice when the satellite passes by the same area about 110 minutes later and

6-15

NDX-120015A

then by thin nilas in the following orbit. If ice concentration is defined strictly as the concentration ofice, assuming a binary surface of either ice and liquid water, even if the lead is 100% covered bygrease ice the estimated ice concentration won't be 100% because the brightness temperature of grease
ice is close to that of open water (O) and the error would be close to 10%. If the lead is covered bynilas, the error would be around 5% if the brightness temperature of nilas is halfway between that ofopen water and thick ice. There are actually advantages in having a retrieval algorithm that provides
such results. Areas of divergence and polynya formations are important to identify (Kwok et al.,2007) and an ice concentration map that provides the location and size of these features provides more
information and would be more useful for scientific research.

2.3.3.3

Melt and Divergent Areas
The Bootstrap Algorithm provides reasonably accurate estimates of ice concentration during

the autumn and winter period when the surface of the ice cover is generally dry. As indicated earlier,the error in ice concentration that may be associated with radiometer noise and uncertainties inemissivity and surface temperature may be as low as 3% during the dry surface period. During latespring and summer, there are known ambiguities. At the onset of melt in spring, the sea ice surface
basically lost its true emissivity because of the transformation of the snow cover into that of a surface
with high emissivity, similar to that of a blackbody. This is because the imaginary part of thedielectric constant of snow becomes relatively high when it is slightly wet (about 3% liquid). Thisphenomenon causes an increase in the error of a few percent because of the higher variability in theemissivity of 100% ice. The one that causes the most serious problem, however, is meltponding
effects. The fraction of meltponded-ice covered surfaces in the Arctic during summer has been
estimated to be as high as 30%. Meltponding is a problem to the algorithm because the surface of
meltponded areas have signatures of ice free water. This makes the tie point for the regions withmeltponding unpredictable and depends on the fraction of meltponds within each data element.
Generally, this causes a bias as much as 30% because of the reduction of the average emissivity of theice cover.

6-16

NDX-120015A

Figure 5.

Summer ice cover in the Arctic as depicted by (a) MODIS data and (b) AMSR-E ice

concentration data.

Color scale for the ice concentration data is the same as in Figure 3.

The problem with meltponding is illustrated in the set of images in the Arctic on 15 June2010 as shown in Figure 5.

The figure to the left is a depiction of the ice cover using the 250m

resolution MODIS 0.65 µm data.

It is apparent in the cloud free portion of the MODIS image thatthe albedo of the area labeled A is much lower (darker) than that in B.

The difference in albedo islikely meltponding which is prevalent in the region during this period.

Such difference in surface

properties is reflected in the ice concentration map for the same day as shown in Figure 5b.

In this

ice map, the A area is shown to have significantly lower concentration than that in the B area.

Inthe images, it is apparent that the area C is likely a meltponded area as well but this is not confirmedin the MODIS image because of cloud cover.
The resolution of the MODIS image in Figure 5 is apparently not good enough to identify
individual melt ponds.

This is in part because the size of individual meltponds is usually too small

to be resolved by the sensor.

Figure 6 shows examples of meltponded areas at a relatively higher

resolution (30m) using Landsat-7 data during an early stage (20 July) and a later stage (8 August)over approximately the same region in the Arctic. These images provide more details about

6-17

NDX-120015A

meltponded areas than the MODIS image but still individual meltponds are not easy to identify.
Aircraft photos over approximately the same time periods are provided in Figure 6 and it is apparent
that the features can be relatively small and sometimes linear.

It is also apparent that different

stages of meltponding can have different overall signature in the passive microwave data.

Figure 6. Meltponding as revealed buy Landsat 7 during early (20 July) and late (8 August) stages.
The aircraft photos were taken at approximately similar periods but for a different year.

A special adjustment of the slope and offset of the dash red line in Figure 1 that determines the
100% line (i.e., AD use for the tie points) through regression analysis as discussed earlier, is done
during the meltponding period but this does not completely solve the problem as illustrated in Figure
5b where the relatively low ice concentration areas do not necessarily represent higher fraction of open
water areas. Comparative analysis of AMSR-E ice concentration data with manually interpreted
MODIS data suggests that the residual error can be as large as 20% in some areas. The date of onsetof melt can be derived from the passive microwave data and the number of days it takes to form
meltponds can be predicted by a thermodynamic model in conjunction with the knowledge of air
temperature. Having this additional knowledge will enable improvements in the determination of

6-18

NDX-120015A

algorithm tie points through more accurate estimates of the percentage of meltponding in eachobservational area.

2.3.3.4

Marginal Ice Zones
Another source of systematic and residual errors is in the characterization of the marginal ice

zones. The general characteristics of the ice edges and adjacent regions was studied by Worby andComiso (2004) using AMSR-E and SSM/I data in conjunction with ship and high resolution satellite
data. Using MODIS data, Comiso and Nishio (2008) show that the ice edge as provided by AMSR-E
is accurate to within about 6 km. With even better resolution than AMSR-E, the AMSR2 ice edge
will likely improve to 5 km. However, there are mismatches between those identified as ice edge by
AMSR-E and those by SSM/I data by as much as 12 km that are associated with the differences inresolution and side lobes. Because of the large contrast in the emissivity of sea ice and open water,
the side lobes cause the brightness temperature of the same surface to be different when the satellite
crosses the ice edge from the ice covered region than when the ice edge is crossed in the opposite
direction. Accurate and consistent identification of the ice edge is important for studies ofinterannual changes and trends in ice extent, especially if data from different microwave sensors are
combined together to form the long-term data record. Such consideration is especially relevant inlight of observed rapid changes in the sea ice extent in the Arctic.
An example of a comparative study of ice cover in the marginal ice zone in the Weddell
Sea using MODIS, AMSR-E and SSM/I data is shown in Figure 7.

The relatively cloud free

MODIS data of ice covered area, as presented in Figure 7a, show relatively high concentration areas(white) inside the ice pack while low concentration values (darker gray values) are found at the ice
edges and in lead and polynya areas. An ice concentration map gridded at 12.5 km resolution derived
from AMSR-E data using the Bootstrap Algorithm is shown in Figure 7b while a similar map using
SSM/I data gridded at 25 km resoluion is presented in Figure 7c. The higher resolution image(Figure 7b) shows some of the mesoscale characteristics of the ice cover but not the spatial details
provided in 7a.

The general features of the ice cover as revealed by AMSR-E and SSM/I data areconsistent except in the marginal ice zone likely due to the differences in resolution and antenna side
lobes.

For quantitative comparisons including MODIS data, plots of data along the transect AB for

the three data sets are presented in Figure 7d.

It is apparent that for this particular location, the ice

distributions as revealed by AMSR-E and SSM/I are very similar with the AMSR-E data showing a
sharper gradient than the SSM/I data. The exact location of the ice edge as determined from theplots for the three sensors appears to be very close to each other and hence a general agreement.
The MODIS data represent surface albedo which fluctuates a lot more than the others in ice covered
areas primarily because of higher sensitivity to changes in surface characteristics.

6-19

NDX-120015A

Figure 7.

Images of the Weddell Sea ice cover in the Southern Hemisphere as depicted by (a)MODIS radiance; (b) AMSR-E ice concentration; and (c) SSM/I ice concentration data.
line is drawn in all 3 images to indicate the location of the 15% ice edge.

A black

The distribution ofradiances from MODIS 0.6 μm channel, AMSR-E ice concentration and SSM/I ice concentrations
along the transect between A and B (in the MODIS) image is shown in (d).

4. Algorithm Implementation
4.1.

Implementation
The discussion of the theoretical basis of the Bootstrap Algorithm in the previous section

provides a general overview of the technique for estimating sea ice concentration and what it takes toobtain accurate values. To summarize, the basic equation for determining sea ice concentration is
given by equation (1) while the associated technique for obtaining the key parameters needed in the

6-20

NDX-120015A

determination is illustrated in Figure 1. The technique is quite simple and should provide accurate
ice concentration values as long as the assumptions for obtaining the “tie-points” for 100% and 0%
ice concentrations are valid.

This is generally the case during dry/cold surface conditions but in thesummer, the uncertainties are increased as discussed earlier.

4.2.

Input/Output parameters:

The input parameters are brightness temperatures as observed by

AMSR2 for each data element.

For ice concentration calculations, these include brightness

temperatures at 18.7 GHz at vertical polarization, 36.5 GHz at both vertical and horizontal
polarization, and 23.8 at vertical polarization.

Input parameters include open ocean tie point,

slope and offset of the regression line (red dash line in Figure 1), and thresholds for oceanmask (see Figure 3).

4.3.

Ancillary data:

The output is ice concentration

Ancillary data includes SST data as supplementary ocean mask derived fromAMSR2 data and land mask derived from published land boundary data supplemented byhigh resolution satellite data.

4.4.

Processing flow:

Although only 1 set of 2 channels are required, 2 sets of channels are usedby the Bootstrap Algorithm to optimize accuracy in the estimate of ice concentration:

VH36

set, representing 36.5 GHz channels at both vertical and horizontal polarization and V1836 set,
representing the 18.7 and 36.5 GHz channels at vertical polarization.

Ice concentration isestimated using the VH36 set for data elements that are above the line AD-4K (see scatter plot
in Figure 1).

The ice concentration for the rest of the data elements are estimated using theV1836 set.

In the process, the data element is also tested to check whether it is in land area

or in the ice free open ocean area.

If it is in either land or ice free ocean, ice concentration is

not calculated and the data is flagged as either land or ocean.

4.4.1

General process to calculate Ice Concentration:

Enter LOOP for each Data Point. Perform 3 tests before attempting to calculate
Ice Concentration:

1 - Test for Invalid Brightness Temperatures
2 - Test for Land
3 - Test if Location is Out-of-Range

If tests call for, set Ice Concentration value and respective Quality Flag(s). If data isfound to be proper, calculate Ice Concentration by hemisphere.
Quality Flags:

1 - a '1' denotes no calculation took place
2 - a '1' denotes invalid brightness temperature
3 - a '1' denotes land location
4 - a '1' denotes latitude is out of ice range
5 - a '1' denotes pixel is out of revised (Dr.Comiso's) land mask area

6-21

NDX-120015A

6 - a '1' SST Filter to cancel the bad IC data
4.4.2

Decision Tree for Ice Concentration:a) Water/Weather Filter Check:watchk1=(weather_slope1 * V22(j,i)) + weather_offset1
watchk2=V22(j,i) - V18(j,i)
watchk3=(weather_slope2 * V37(j,i)) + weather_offset2
watchk4=(VH37_threshold_slope * V37(j,i)) + VH37_threshold_offset
if ((watchk1 .gt. V18(j,i) .or. watchk2 .gt. weather_threshold .or.
(watchk3 .gt. V06(j,i) .and. V06(j,i) .gt. mintb) .and.
(H36(j,i) .le. watchk4 .or. V36(j,i) .ge. 230.0)) then ICE_DATA(j,i)=0.0b) Channel Pair Check:vh37chk=(VH37_100%_iceline_offset-8%offset)+(VH37_100%_iceline_slope * V36(j,i))if (H36(j,i) .gt. vh37chk) then"Calculate ice using V36 & H36 channels"
else"Calculate ice using V36 & V18 channels"
endif

4.4.3

Parameters for Ice Concentration Algorithm:

Sample routine to retrieve Arctic input parameters:
SUBROUTINE ICE_PARA_SUB_N(month,day,wtp,itp,>

wintrc,wslope,wxlimt,wintrc2,wslope2,>

vh37,v1937)

INTEGER*2 month,day,season
REAL*4 wtp(3),itp(3)REAL*4 wintrc,wslope,wxlimt,wintrc2,wslope2
REAL*4 vh37(2),v1937(2)SEASON=1
IF (MONTH .GE. 6 .AND. MONTH .LE. 9) SEASON=2
IF (MONTH .EQ. 10 .AND. DAY .LE. 15) SEASON=2
IF (SEASON.EQ.1) THEN
WINTRC= 83.73
WSLOPE= 0.5352
WXLIMT= 18.39

6-22

NDX-120015A

WINTRC2= 12.22
WSLOPE2= 0.7020
ELSE
WINTRC= 82.71
WSLOPE= 0.5352
WXLIMT= 23.34
WINTRC2= 12.22
WSLOPE2= 0.7020
ENDIF
WTP(1)= 207.2
WTP(2)= 131.9
WTP(3)= 182.4
ITP(1)= 256.3
ITP(2)= 241.2
ITP(3)= 258.9
IF (VH37(1) .eq. -999) THEN
VH37(1)= -30.26
VH37(2)= 1.0564
ENDIF
IF (V1937(1) .eq. -999) THENV1937(1)= 110.03
V1937(2)= 0.5759
ENDIF
RETURN
END

Sample routine to retrieve Antarctic input parameters:
SUBROUTINE ICE_PARA_SUB_S(month,day,wtp,itp,>

wintrc,wslope,wxlimt,wintrc2,wslope2,>

vh37,v1937)
INTEGER*2 month,day
REAL*4 wtp(3),itp(3)REAL*4 wintrc,wslope,wxlimt,wintrc2,wslope2
REAL*4 vh37(2),v1937(2)WINTRC= 85.13
WSLOPE= 0.5379

6-23

NDX-120015A

WXLIMT= 18.596
WINTRC2= 10.93
WSLOPE2= 0.7046
WTP(1)= 207.6
WTP(2)= 131.9
WTP(3)= 182.7
ITP(1)= 259.4
ITP(2)= 247.3
ITP(3)= 261.6
IF (VH37(1) .eq. -999) THEN
VH37(1)= -38.31
VH37(2)= 1.0969
ENDIF
IF (V1937(1) .eq. -999) THENV1937(1)= 114.26
V1937(2)= 0.5817
ENDIF
RETURN
END

4.5.

Examples of output using AMSR-E data are presented in Figures 4, 5 and 7.

See discussion

in the previous section for details.

4.6.

The basic limitations of ice concentration data from AMSR2 include relatively coarse
resolution and the inability to unambiguously identify areas covered by thin ice, pancake ice
and meltponded ice.

Thin ice, pancake ice and meltponded ice can have emissivities

intermediate to those of open water and dry thick ice and could contribute to significant
uncertainties in the retrieval of ice concentration.

A good scheme that enables classificationof each data element into different surface types would help minimize uncertainties.

5. Validation Concept
Validation of ice concentrations from AMSR2 will be done primarily using high-resolution
satellite data as has been done previously (Comiso and Steffen, 2001).

Aircraft data from the Ice

Bridge Project over both Northern and Southern Hemispheres will also be utilized for directcomparison and to validate interpretation of high-resolution satellite data.

The data include detailed

topography information and high resolution photos of the surface twice a year over key sea ice
regions of the Arctic and the Antarctic.

Special emphasis will be done on retrieved ice

6-24

NDX-120015A

concentrations over new ice, marginal ice zones, and meltponded regions.

References
Bhatt, U.S., D. A. Walker, M. Raynolds, J. C. Comiso, et al., 2010: Circumpolar Arctic tundra
vegetation change is linked to sea-ice decline, Earth Interactions, 14, 1-20,
doi:10.1175/2010EI315.1.
Broecker, W. S. 1991: The Great Ocean Conveyor, Oceanography, 4, 79-89.
Cavalieri, D.J., P. Gloersen, W.J. Campbell, 1984: Determination of sea ice parameters with the
Nimbus7 SMMR, J. Geophys Res., 89, 5355-5369, 1984.
Cho, K., N. Sasaki, H. Shimoda, T. Sakata and F. Nishio, 1996: Evaluation and improvement of
SSM/I sea ice concentration algorithms for the Sea of Okhotsk, J. Remote Sensing ofJapan, 16(2), 47-58, 1996.
Colony, R. and A. Thorndike, 1985: Sea ice motion as a drunkard’s walk. J. Geophys. Res., 90965-974.
Comiso, J. C., 1983: Sea ice microwave emissivities from satellite passive microwave andinfrared observations, J. Geophys. Res., 88(C12), 7686-7704.
Comiso, J. C., 1986: Characteristics of winter sea ice from satellite multispectral microwave
observations, J. Geophys. Rev., 91(C1), 975-994, 1986.
Comiso, J.C., 1985: SSM/I Concentrations using the Bootstrap Algorithm, NASA RP,
1380, 40pp.
Comiso, J. C., 2002: A rapidly declining Arctic perennial ice cover, Geophys Res. Letts.,
29(20), 1956, doi:10.1029/2002GL015650.
Comiso, J. C., 2009:

Enhanced sea ice concentration and ice extent from AMSR-E, J. Remote

Sensing Soc. of Japan, 29(1), 199-215.
Comiso, J. C., (2010), Polar Oceans from Space, Springer, New York, 518pp.
Comiso, J. C. and F. Nishio, 2008: Trends in the sea ice cover using enhanced and compatible
AMSR-E, SSM/I, and SMMR data, J. Geophys. Res. 113, C02S07,
doi:10.1029/2007JC004257.
Comiso, J. C. and C. L. Parkinson, 2008: Arctic sea ice parameters from AMSR-E using two
techniques, and comparisons with sea ice from SSM/I, J. Geophys. Res., 113, C02S05,
doi:10.1029/2007JC004255.
Comiso, J.C., and K. Steffen, 2001: Studies of Antarctic sea ice concentrations from satellite
data and their applications, J. Geophys. Res., 106(C12), 31361-31385.
Comiso, J. C., D. J. Cavalieri, and T. Markus, 2003. Sea ice concentration, ice temperature, andsnow depth, using AMSR-E data, IEEE TGRS, 41(2), 243-252.
Comiso, J.C., C.L. Parkinson, R. Gersten, and L. Stock, 2008: Accelerated decline in the Arctic

6-25

NDX-120015A

sea ice cover, Geophy. Res. Lett. 35, L01703, doi:10.1029/2007GL031972.
Comiso, J.C., D. Cavalieri, C. Parkinson, and P. Gloersen, 1997: Passive microwave algorithms
for sea ice concentrations, Remote Sensing of the Env., 60(3), 357-384.
Eicken, H., M.A. Lange, and G.S. Dieckmann, 1991: Spatial variability of sea-ice properties in
the northwestern Weddell Sea, J. Geophys. Res., 96, 10,603-10,615.
Gloersen P., W. Campbell, D. Cavalieri, J. Comiso, C. Parkinson, H.J. Zwally, 1992: Arctic andAntarctic Sea Ice, 1978-1987: Satellite Passive Microwave Observations andAnalysis, NASA Spec. Publ. 511.
Grenfell, T.C. 1992.Surface-based passive microwave studies of multiyear ice. J.
Geophys. Res., 97(C3), 3485-3501.
Kumerow,

C, 1993: On the accuracy of the Eddington approximation for radiative transfer inthe microwave frequencies," J. Geophys. Res, Vol. 98, pp. 2757- 2765.
Kwok, R., J. C. Comiso, S. Martin, and R. Drucker, 2007: Ross Sea Polynyas: Response of ice
concentration retrievals to large areas of thin ice, J. Geophys. Res., 112, C12012,doi:10.1029/
2006JC003967.
Lindsay, R.W., and J. Zhang, 2005:

The thinning of Arctic sea ice , 1988-2003: Have we reached

the tipping point?, J. Clim., 18, 4879-4894.
Markus, T and D. Cavalieri, 2000: An enhancement of the NASA team sea ice algorithm, IEEE
Trans. Geosci. Remote Sensing, 38, 1387-1398.
Martin, S., R. Drucker, R. Kwok, and B. Holt, 2005: Estimation of the thin ice thickness and heat
flux for the Chukchi Sea Alaskan coast polynya from Special Sensor Microwave Imager data,
1990-2001, J. Geophys. Res, 109, C10012, doi:10.1029/2004JC002428.
Matzler, C., R. O. Ramseier, and E. Svendsen, 1984: "Polarization effects in sea ice
signatures," IEEE J. Oceanic Engineering, Vol. OE-9, pp. 333-338, 1984.
Naoki, K., J. Ukita, F. Nishio, M. Nakayama, J.C. Comiso, and A. Gasiewski (2008), Thin sea ice
thickness as inferred from passive microwave and in situ observations, J. Geophys. Res.
113, C02S16,doi:10.1029/2007JC004270.
Parkinson, C. L., and J. C. Comiso, 2008:, Antarctic sea ice from AMSR-E from two algorithms
and comparisons with sea ice from SSM/I, J. Geophys. Res., 113, C02S06,
doi:10.1029/2007JC004253.
Parkinson, C. L., J. C. Comiso, H. J. Zwally, D. J. Cavalieri, P. Gloersen, and W. J. Campbell,
1987: Arctic Sea Ice 1973-1976 from Satellite Passive Microwave Observations, NASA Spec.
Publ. 489.
Serreze, M.C., M.M. Holland, and J. Stroeve, 2007: Perspective on the Arctic's shrinking sea-ice
cover, Science, 315, 1533-1536.
Shibata, A., H. Murakami, and J. Comiso, 2010: Anomalous Warming in the Arctic Ocean in the

6-26

NDX-120015A

Summer of 2007, J. Remote Sensing Society of Japan, 30(2), 105-113.
Smith, W.O. Jr., and J. C. Comiso, 2009: Southern ocean primary productivity: Variability and aview to the future, in Smithsonian at the Poles: Contributions to the International Polar Year
Science, eds. I. Krupnik, M.A. Lang and S.E. Miller, Smithsonia Institution Scholarly Press,
Washington, DC.,

pp. 309-318.

Steffen, K., D. J. Cavalieri, J. C. Comiso, K. St. Germain, P. Gloersen, J. Key, and I.
Rubinstein, 1992: "The estimation of geophysical parameters using Passive Microwave
Algorithms," Chapter 10, Microwave Remote Sensing of Sea Ice, (ed. by Frank
Carsey), American Geophysical Union, Washington, D.C., 201-231.
Svendsen, E., C. Matzler, T.C. Grenfell, 1987: A model for retrieving total sea ice
concentration from a spaceborne dual-plarized passive microwave instrument
operating near 90 GHz, Int. J. Rem. Sens., Vol. 8, pp. 1479-1487.
Swift, C.T., L.S. Fedor, and R.O. Ramseier, 1985: An algorithm to measure sea ice
concentration with microwave radiometers, J. Geophys. Res., 90(C1), 1087-1099.
Tucker, W.B., D.K. Perovich, and A.J. Gow, 1992: "Physical properties of sea ice relevant toremote sensing," Chapter 2, Microwave Remote Sensing of Sea Ice, (ed. by Frank
Carsey), American Geophysical Union, Washington, D.C., 9-28.

Turner, J., J.C. Comiso, G. J. Marshall, T.A. Lachlan-Cope, T. Bracegirdle, T. Maksym,
M. Meredith and Z. Wang, 2009: Non-annular atmospheric circulation change induced bystratospheric ozone depletion and its role in the recent increase of Antarctic sea ice extent,
Geophy. Res. Lett. 36, L08502, doi:10.1029/2009GL037524
Vant, M.R., R.B. Gray, R.O. Ramseier, and V. Makios, 1974: Dielectric properties offresh and sea ice at 10 and 35 GHz, J. Applied Physics, 45(11), 4712-4717.
Wang, M. and J. Overland 2009: A sea ice free summer Arctic within 30 years?, Geophys. Res.
Lett., 36, L07502, doi:10.1029/2009GL037820.
Wadhams, P. 2000: Ice in the Ocean, Gordon and Breach Science Publishers, London.
Weeks, W.F., and S. F. Ackley, 1986: The growth, structure and properties of sea ice, The
Geophysics of Sea Ice, edited by N. Unterstiener, pp. 9-164, NATO ASI Ser.B, vol.
146, Plenum, New York.
Worby, A. P., and J. C. Comiso, 2004: Studies of Antarctic sea ice edge and ice extent fromsatellite and ship observations, Remote Sensing of the Environment, 92(1), 98-111.
Yaguchi, R. and K. Cho, 2009: Validation of sea ice drift vector extraction from AMSR-E and SSM/I
data using MODIS data, J. Remote Sens. Japan, 29(1), 242-252.
Zwally, H.J., J.C. Comiso, C. Parkinson, D. Cavalieri, P. Gloersen, 2002: Variability of the
Antarctic sea ice cover, J. Geophys. Res. 107(C5), 1029-1047.
Zwally, H. J., J. C. Comiso, C. L. Parkinson, W. J. Campbell, F. D. Carsey, and P.

6-27

NDX-120015A

Gloersen, 1983: Antarctic Sea Ice 1973-1976 from Satellite Passive Microwave
Observations, NASA Spec. Publ. 459.

6-28

NDX-120015A

Chapter 7.
Description of GCOM-W1 AMSR2Snow Depth Algorithm

Richard Kelly
Department of Geography & Environmental Management
University of Waterloo
7-1

NDX-120015A

Table of Content

1.

Introduction ................................................................................................................... 3

2.

Algorithm Overview ....................................................................................................... 3

3.

Theoretical Description .................................................................................................. 4

4.

Algorithm Implementation ............................................................................................. 6

5.

4.1.

Implementation ......................................................................................................... 6

4.2.

Input/Output parameters ............................................................................................ 8

4.3.

Ancillary data ........................................................................................................... 8

4.4.

Processing flow ........................................................................................................ 9

4.5.

Example output......................................................................................................... 9

4.6.

Limitations ............................................................................................................. 10

Validation Concept ....................................................................................................... 13

References ........................................................................................................................... 14

7-2

NDX-120015A

1. Introduction
Seasonal snow cover is a dynamic hydrological variable in the Earth’s water cycle with changes
in northern hemisphere snow cover extent ranging from a maximum of about 47 x 106 km2(including the Greenland Ice Sheet) in January and February to 4 x 106 km2 in August (Robinson andFrei, 2000). From a thermodynamics perspective, this change in land surface albedo and subsequent
can have important implications for atmospheric processes with changes in Siberian high pressure
system tracks having been attributed to the early onset of a seasonal snow cover in Eurasia (Gong etal., 2006). Assessing global snow accumulation is arguably a more challenging task in themonitoring of seasonal snow dynamics. The use of in situ measurements to estimate snow waterequivalent or snow depth is subject to complex scaling-up assumptions are often uncertain,
especially when the measurements are sparsely distributed. Satellite-based estimates of snow
accumulation (snow depth, or SD, and snow water equivalent, or SWE) using microwave
observations are offer significant potential since measurements are consistent, frequent (twice daily
in many high latitude regions), are capable of day and night observations and are able to observe theEarth’s surface under all but the thickest cloud conditions. Several studies have shown that passive
microwave observations of snow can yield important information about snow climatology but havealso been subject to large uncertainties in SD or SWE estimates from the satellite instruments used
(e.g. Biancamaria et al., 2008; Brown et al., 2010; Wulder et al., 2007). Recent advances in retrieval
science for passive microwave estimates of snow have been reported in the literature (e.g. Kelly etal., 2009). Recently, Derksen (2008) demonstrated the important contribution that the low frequency
passive microwave channels can make to SWE estimates. Such studies are important and help toreduce uncertainties in the practical application of microwave radiometry to estimate SWE and SD.
The implementation of a global algorithm to estimate SD constitutes a significant challenge. This
document describes the approach taken to estimate SD using the Japanese Aerospace Exploration
Agency (JAXA) Advanced Microwave Scanning Radiometer 2 (AMSR2) which is the main payload
aboard the Global Change Observation Mission 1st Water (GCOM-W1). This document, an
Algorithm Theoretical Basis Development (ATBD) document provides a description of thealgorithm science as implemented for AMSR2. It is also a follow on from the AMSR-E approach.

2. Algorithm Overview
The GCOM-W1 AMSR2 instrument constitutes a key component in the record of satellitepassive microwave observations of Earth from space. Continuing on from the Advanced Microwave
Scanning Radiometer for EOS (AMSR-E) launched in 2002, the Special Sensor Microwave Imager(SSM/I) launched in 1987 and the Scanning Multichannel Microwave Radiometer (SMMR)
launched in 1979, AMSR2 will extend the measurement record beyond 35 years enabling significant
7-3

NDX-120015A

insight into snow accumulation variability at regional to global scales. The current algorithm
approach, therefore, has a core aim to leverage maximum information and knowledge from the
AMSR2 observations and provide insight into the nature of snow accumulation changes at regional
to global scales over the AMSR2 lifetime and back through the AMSR-E, SSMI and SMMR
observation eras.
The current AMSR2 SD algorithm is an evolution of the original AMSR-E SWE algorithm andtakes advantage of the expanded range of channels available on the AMSR-E instrument comparedwith the SSM/I and SMMR. In principle, the current algorithm is a frequency difference algorithm
approach that builds on the work of Chang et al. (1996). It undertakes a forest correction, retrieves
shallow snow using the 89 GHz channel, and retrieves moderate snow accumulation using the 37
GHz channel. In addition, a deep snow estimation method is implemented using the 19 and 10 GHz
channel. To exploit the native resolution of the AMSR2 instrument, retrievals are conducted at theoriginal native spatial resolutions of the channels and then gridded to 0.25x0.25 degree and
0.10x0.10 degree grids. Using the study of Dewey and Heim (1986) locations where snow
accumulation is climatologically unlikely are masked out (e.g. at low latitudes and over the oceans)
and retrievals are not conducted over permanent terrestrial ice surfaces (ice sheets and glaciers). The
algorithm is calibrated for SD because there are many more SD measuring sites in the northernhemisphere for testing the depth estimates than for testing SWE.

3. Theoretical Description
The upwelling radiation measured at a satellite passive microwave radiometer is the apparent
temperature (Tb) and can be expressed as:
𝑇𝑇 = �𝑅 𝑇𝑇𝑠𝑠𝑠 + (1 − 𝑅)𝑇𝑇𝑠𝑠𝑠𝑠 �𝑒 −𝜏 + 𝑇𝑇𝑎𝑎𝑎

where e-τ is the atmospheric transmissivity R is the surface reflectivity,

Tb sky(1)

is the sky radiation,

Tb surf is the surface emission and Tb atm is the atmospheric (Chang et al., 1996). In general Tb atm andTb sky are very small and are ignored. Therefore, Tb is directly related to surface features and if
microwave atmospheric windows are used, atmospheric transmissivity is maximized (~1).
The wavelength-dependent intensity of radiation emitted from a snowpack depends on several
physical variables: kinetic temperature of the snow, snow grain size distribution and the grainvolume to air fraction, the underlying surface conditions and in situ vegetation characteristics. At thesatellite scale, the fraction of water within the instantaneous field of view is also important (Gunn etal., 2011). For dry snowpacks with thickness greater than 5 cm or SWE greater than 10 mm,
isotropic scattering of naturally emitted microwave radiation by snow crystals occurs controlling the
7-4

NDX-120015A

emission from a dry snowpack (Ulaby and Stiles, 1980). Below the 20 GHz frequency threshold,snow is almost transparent. For remote sensing of snow applications, this emission behavior is
detected at frequencies greater than about 20 GHz (e.g., Chang et al., 1987) and has formed thebaseline approach for SD and SWE retrievals.
Satellite passive microwave instruments have been measuring natural upwelling microwaves
from the Earth for 33 years. Generally, for snow covered terrain the strength of scattering signal is
proportional to the SWE or snow depth, and it is this relationship that can be used for estimating thewater equivalent. Several algorithms for the retrieval of SWE or SD have been proposed in theliterature, such as those by Aschbacher (1989), Chang et al. (1987) reviewed and updated forforested areas by Foster et al. (1997). These approaches typically use the brightness temperature(Tb) difference between 19 and 37 GHz as an estimate of scattering strength which is the calibrated
statically to SWE. More recent studies have shown an improvement to the estimation of SWE atregional scales through the use of numerical techniques for the inversion of semi-empirical
relationships (e.g. Tedesco et al., 2004). Grippa et al. (2004) have also used a dynamic algorithm toestimate seasonal SWE in Siberia, although this approach is more for climatological rather thanhydrology application studies. A dynamic algorithm developed by Kelly et al. (2003) demonstrated
that the original statically-parameterized algorithms (Chang et al., 1987 and Foster et al., 1997) canbe improved upon by explicitly accounting for snowpack changes in physical state (e.g. grain sizeand density). The current standard JAXA SD algorithm implemented for AMSR-E has evolved fromthis approach but with parameterization of grain size based on Tb polarization signatures (Kelly,2009). This is also the approach that forms the baseline implementation for the AMSR2
implementation.
Within the regions of snow possibility, snow is tested for as either moderate to deep snow or as
shallow snow. For moderate to deep snow to be present, there are two criteria that have to be met.
First, the difference between non-scattering and scattering channel must be greater than zero.
Historically the following equation is often used:
𝑇𝑇18𝑉 − 𝑇𝑇36𝑉 > 0(2)

where Tb 18V and Tb 36V are the vertically-polarized brightness temperature at 18 GHz and 36 GHz
respectively. This approach has demonstrated that a positive brightness temperature difference
between a ‘non-scattering channel (18 GHz) and a scattering channel (36 GHz) indicates thepresence of a scattering medium (Chang et al., 1987). Scattering effects can also be observed at the
18 or 19 GHz wavelengths for deeper snowpacks (Tsang et al., 2000). The use of the 10 GHz
channel on AMSR2 provides greater insensitivity to SD. For moderate to deep snow accumulations
to be detected, either of the following two conditions must be satisfied:

7-5

NDX-120015A

𝑇𝑇10𝑉 − 𝑇𝑇36𝑉 > 0(3)

𝑇𝑇10𝐻 − 𝑇𝑇36𝐻 > 0(4)

where Tb 10V and Tb 10H are the vertically-polarized and horizontally-polarized brightness
temperatures at 10 GHz respectively, and Tb 36V and Tb 36H the same but for 36 GHz.
The dominant physical mechanism used for estimating the presence and the amount of snow isvolume scattering. Scattering can be caused not only by dry snow but also by rainfall hydrometeors
and so a threshold must be identified to separate out these two scattering types (e.g. Grody andBasist, 1996). Typically, rainfall has a higher Tb at 36 GHz than snow suggesting that the rainfall
distribution in the feature space plots should be located at higher Tb 36V values greater than snow
distribution. At Tb 36V , the threshold is 255K whilst at Tb 36H it is 245K (Kelly, 2009)
A shallow snow detection capability for AMSR2 is also applied in the same way as that for
AMSR-E. Shallow snow is particularly challenging to detect, since it is mostly transparent at Tb 36
frequencies. By using the combination of the 23 and 89 GHz channels at vertical and horizontal
polarizations, and a surface temperature estimator, T phys, devised by Kelly et al. (2003):
𝑇𝑝ℎ𝑦𝑦 = 58.08 − 0.39𝑇𝑇18𝑉 + 1.21𝑇𝑇23𝑉 − 0.37𝑇𝑇36𝐻 + 0.36𝑇𝑇89𝑉

(K)

(5)

T phys, is accurate to ±6 K. Locations where T phys, is less than 267 K are flagged as locations whereshallow snow is possible if moderate to deep snow has not been previously detected.

To ensure thatshallow snow is present, the following condition is evaluated:
𝑇𝑇23𝑉 > 𝑇𝑇89𝑉 & 𝑇𝑇89𝑉 < 255𝐾 & 𝑇𝑇23𝐻 > 𝑇𝑇89𝐻 & 𝑇𝑇89𝐻 < 255𝐾

(6)

If the condition is true then shallow snow is present. The 255K threshold is a conservative value toensure that only cold Tbs are used that more probably relate to snow cover emission. By ensuringthat Tb 89V and Tb 89H are less than the 255 K threshold in (6) and that the surface temperature (5) isless than 267 K the probability that atmospheric contamination is present is minimized. In this way,
the 89 GHz measurements are used conservatively.

4. Algorithm Implementation
4.1.

Implementation

The current standard retrieval is performed using L1R Global Swath data . SD retrievals are
performed on the instantaneous field of view (IFOV) samples. The procedure followed is described
7-6

NDX-120015A

below:
Step 1. Obtain Tbs. AMSR2 brightness temperatures are acquired from L1R data and used at native
channel resolution with the exception of the 89 GHz channel, which is re-sampled to the 37 GHz
channel resolution.
Step 2. Calculate land surface physical temperature. The surface physical temperature (T) in
Kelvin is from Kelly et al. (2003) where Tb 18V , Tb 22V , Tb 36H and Tb 89V are Tb FFP brightness
temperatures at FF frequency (GHz) and P polarization (Vertical or Horizontal):

T = 58.08 - 0.39Tb 18V + 1.21Tb 22V - 0.37Tb 36H + 0.36Tb 89V

[K]

(7)

Step 3. Test for moderate to deep snow presence. Thresholds are checked to ensure cold snow
conditions are potentially present in the 36 GHz brightness temperatures (Tb 36H <245K &Tb 36V <255K). If the condition is true then snow is possible and shallow or medium depth of snow isretrieved (see Step 5). If this condition is not met then a shallow snow test is performed (step 4).
Step 4. Test for shallow SD. If Tb 10V -Tb 36V > 0 K or Tb 10H -Tb 36H > 0 K then medium to deep snow
is assumed to be present (go to full retrieval, Step 5). Otherwise snow presence is possible but it islikely to be shallow snow if: Tb 89V <= 255 K and Tb 89H <= 265 K & Tb 23V -Tb 89V >0 K &Tb 23H -Tb 89H >0 K & T < 267 K.

If shallow snow is detected, SD is estimated as 5.0 cm.

Step 5. Moderate to deep SD retrieval. For the retrieval of the SD the following general procedure
is followed:

SD = ff(SD f ) + (1-ff)*(SD o ) [cm](8)where SD is the total sample snow depth, SD f is the snow depth from the forest component of theinstantaneous field of view (IFOV) and SD o is the snow depth from non-forested component of the
IFOV:
SD f = 1/log 10 (pol 36 ) * (Tb 18V -Tb 36V )/(1-fd * 0.6)[cm]

SD o = [1/log 10 (pol 36 ) * (Tb 10V -Tb 36V )] + [1/log 10 (pol 18 ) * (Tb 10v -Tb 18V )]

(9)
[cm]

(10)

where ff is the forest fraction (where a ff of 1.0 = 100% forest fraction and a ff of 0.0 = 0% forest
fraction) from the MOD12Q1IGBP product, and fd is the forest high spatial resolution (500m) forest
density from University of Maryland Vegetation Continuous Field (VCF) data circular smoothed at
15km diameter and re-gridded to global 1 km. The polarization factors, pol 36 = Tb 36V -Tb 36H andpol 18 = Tb 18V -Tb 18H are constrained through optimization to vary conservatively through the season
as the mean depth-integrated snowpack grain size increases (Kelly, 2009). Note the use of thedifference between the 18 and 36 GHz channels to maximize spatial resolution in forested areas, and

7-7

NDX-120015A

the use of 10V-36V (increased dynamic range) and Tb 10V -Tb 18V for deep snow. Note also that fd isscaled (0.6) through optimization of validation data.

4.2.

Input/Output parameters

The input parameters are divided into two categories: dynamic swath granule brightness
temperature file (L1R) and static ancillary data files designed to characterize forest cover, land,
oceans, coasts and ice cover and snow climatology. JAXA AMSR2 L1R swath data form the basis of
the SD retrievals and are obtained at EORC as part of the processing chain. The output file for the
first step is a swath product of four raster grid arrays. The first two are decimal latitude andlongitude grids (243 x 2100 floating point values), the third array contains snow depth estimates incentimetres (243 x 2100 cell array of floating point values) and the fourth data layer is a cell array offlags that identifies different surfaces or bad data. The following flags are included:
Snow possible:

0

Water:

16

Climatologically snow impossible:

32

Permanent ice:48

Tb out of range:

192

Bad spacecraft attitude:208

Bad Tb:224

The ancillary data sets are outlined in the next section.
The initial algorithm processing step, where the primary retrieval is performed, provides snow
depth values for all samples in the swath granule. Retrievals are conducted on the day and nighttime
data. The retrieval data are retained as swath data (scene data) before the processing suite computedaily global (latitude and longitude) snow depth and north polar stereographical projected snow
depth for day and night. Finally, the data are composited to produce monthly global and polar
projected snow depth estimates.

4.3.

Ancillary data

The following data sets are used to parameterize the retrieval algorithm. These are static data sets
and are created for the algorithm and ingested at algorithm runtime.
•

Global forest fraction from Boston University data (MOD12Q1IGBP) circular smoothed (15 km
diameter) and mapped to 0.00833° grid (Hansen et al., 2003).

•

Global forest ‘spatial density’ from University of Maryland Vegetation Continuous Fraction(VCF) (500m) circular smoothed (15 km diameter) and re-gridded to 0.00833°.

•

Land, Ocean Coasts & Ice mask derived from MODIS MOD12Q1 IGBP land cover data(collection V004). If more than 50 % water cover is present in a 25 x 25 km EASE-Grid cell, the

7-8

NDX-120015A

grid cell is flagged as water. Also, snow climatology data set (Dewey and Heim, 1982) is used todetermine the possibility/ impossibility of snow presence.

4.4.

Processing flow

AMSR2
L1R Tbs

Deg. Lat.

Forest Fraction
MOD12

Deg.
Lon.
Tests for Tb
quality andland covertypes

Forest ‘density’
UMD

Flags

Snow possible

Snow
Depth

AMSR2 Snow
Depth granule(~15 per day)Land, ocean, ice
snow climatology

Snow
depth
retrieval

Daily snow depth
composite: polar

Month snow depth
composite: polar

Daily snow depth
composite: global

Daily snow depth
composite: global

Product data
Figure 1 Processing chain for the AMSR2 snow depth algorithm production.

4.5.

Example output

7-9

An
cill
ar
y
da
ta

NDX-120015A

Figure 2 shows five example outputs extracted from browse data at http://gcom-w1.jaxa.jp web sitefrom where data can be obtained. Each of the five figures shows the output from the algorithm
implementation. Figure 2a shows the granule data and its swath coverage for a pass on 30 March
2011. Figure 2b and 2c show a daily composite of granule data for the same day gridded into a polarstereographic projection and global latitude/longitude projection respectively. Figure 2d and 2e are
the same but composited for the month of March 2011.
a.

c.

b.

d.

e.

Figure 2. Snow depth products for a granule over the western Pacific on 30 march 2011 (a), thecomposited daily granules for 30 March 2011 projected to a polar stereographic (b) and Platte Carre(c) projection, and the composited monthly snow depth for March 2011 in polar stereographic (d)and Platte Carre (e) projections.

4.6.

Limitations

7-10

NDX-120015A

The simplicity of the algorithm implementation allows a fast inversion. However, it is also thesource of errors that are known to be related to the retrieval of snow depth from satellite passive
microwave observations: (i) attenuation of snow emission by atmospheric effects caused by variationsin atmospheric thickness; (ii) attenuation of snow emission by forest and tall stand vegetation; (iii)
evolution of snow grain size and density; (iv) fractional emission effects on mixing by water bodies (infrozen and liquid states); (v) inability to map the snow depth of wet snow.
Atmospheric Attenuation: the microwave signal emitted at the surface passes through theatmosphere before being detected by a space-borne sensor, and thus it is subject to the effects ofatmospheric absorption and emission. This attenuation varies with optical depth of the atmosphere(as a function of variable land surface elevation) and air mass type and microwave frequency. As aconsequence the difference between the 19 and 37 GHz brightness temperatures from space-borne
data can be different from those obtained using brightness temperatures at the surface (Wang andTedesco, 2007; Markus et al. 2006). The Tbs observed by the AMSR-E depend on the effective
atmospheric temperature (T a ), optical depth (τ), surface emissivity (e), and surface temperature (T s ).
Rawinsonde data from different stations can be used to stratify T a and τ. The surface emissivity e canbe derived through an atmospheric absorption model (e.g. Rosenkranz, 1998) and the surface
brightness temperature can then obtained by the product of e and the land surface temperature T s .
This is for clear-sky conditions and is the rationale for the work presented by Savoie et al. (2009)
who simplified this background science and corrected retrievals over the Tibetan plateau by simple
atmospheric Tb corrections based on an empirical model and elevation information. In case ofcloudy skies it is necessary to have information about the cloud thickness and liquid water content.
Recent results show that even under a clear sky the atmospheric absorption could account for asmuch as ∼ 25-50 % to the estimation of snow depth or SWE, depending on the reference brightness
temperatures (Wang and Tedesco, 2007).
Forest Attenuation: Forest cover represents a significant source of error for satellite passive
microwave snow depth and SWE retrieval algorithms. The presence of forest attenuates the radiation
emitted by the underlying snowpack, hence affecting the retrieval accuracy of the algorithm. The
problem is complex: stem volume and canopy closure and gap fraction within a footprint areimportant modulators of the passive microwave emission. At the northern edge of the boreal forest,
stunted conifers that are sparse might be considered forest and yield a reasonably accurate snow
depth estimate. Yet in the middle of the taiga, a pixel covered by dense forest growth would yieldmuch less SWE by passive microwave measurements. With newer high spatial resolution vegetation
products available, such as the Vegetation Continuous Fraction (VCF) product (Hansen et al., 2003)
and the global Landsat GeoCover data base at the University of Maryland, significant improvement
in forest (and other vegetation type) parameterization is possible. Recent work by Metsamaki et al.
(2005) demonstrates and effective approach to mapping snow in forested lands which could be of

7-11

NDX-120015A

benefit to improved forest parameterization.
Snow microphysical properties: Kelly and Chang (2003) demonstrated that SWE retrieval
errors can be caused by not accounting for snowpack physical properties (mostly snow grain size anddensity) which metamorphose during the snow season. Snow crystals metamorphose in response tovapour gradients within the snowpack and as a result of melting and refreezing cycles which alter the
structure of the bulk snowpack. Also snow bulk density changes during the season, usually
increasing. These physical properties control the microwave response from the snow (Armstrong etal., 1993). Although empirical and physically-based models have been developed to predict thegrowth of the snow crystals and density (e.g. Jordan, 1991) it is not straightforward to select ageneral model that will account for regional to global scale conditions.
The current ASMR-E depth algorithm uses polarization information to attempt to account for the
temporal evolution of a snowpack (especially grain size). The polarization difference at 36 GHz (and19 GHz for thicker snowpacks) changes as a snowpack grain size metamorphoses, as demonstrated
in field studies by Matzler (1987). However, it is known that the current standard algorithm’s
performance is substandard in the north east Siberian regions as demonstrated by Clifford (2010) andis likely caused by the grain size parameterization in (3) and (4) requiring improvement.
Effect of terrestrial water bodies: Lakes and ponds affect the accuracy of passive microwave
snow depth and SWE retrievals because the strong gradient in brightness temperatures between
liquid and frozen water changes the Tb difference relationship between 19 and 37 GHz channels.
Analysis of airborne passive microwave data acquired in the Northwest Territories, Canada in April
2005 showed the relationship between Tb 37V and lake cover fraction is reversed across the northern
boreal forest compared to the open tundra (Derksen et al., 2005).

Over forested sites, lower Tb 37V

measurements were observed over lakes relative to land, while the Tb 19V data showed little
sensitivity to lakes. Conversely, at tundra sites the Tb 37V values were greater over lakes than overterrestrial surfaces. This difference in response to lake ice at Tb 37V impacts SWE retrievals because
the increase in brightness temperature at 37 GHz across lake rich tundra areas reduces Tb 37V – Tb 19V
differences. Tacit acknowledgement of this process by Derksen et al. (2010) led to a tundra specificalgorithm for SWE estimation based on a time series approach.
Wet snow mapping: The presence of liquid water within the snowpack increases the absorption
and emission of radiation, reducing the emission depth from the surface. As a consequence the Tb
difference at 19 and 37 GHz approaches 0 K completely terminating the sensitivity to depth or SWE.
It is important, therefore, to distinguish between dry and wet snow conditions, to exclude thosepixels containing wet snow in order to reduce the uncertainty on the SWE retrieval. Walker andGoodison (1993) reported a wet snow discrimination technique based on the 37 GHz polarizationdifferences while Sun et al. (1996) used a neural network to determine snow wetness in vegetatedterrain. Identification of wet snow, therefore, is an important aspect of snow depth and SWE

7-12

NDX-120015A

mapping with passive microwave observations.

5. Validation Concept
Figure 3 shows the locations of 242 and 254 snow depth recording sites during the 2002-2003
and 2003-2004 winter seasons that are from the World Meteorological Organization’s (WMO) snow
depth data archive. These data were compared with AMSR-E snow depth estimates from the current
algorithm and the original Chang et al. (1987) algorithm. Table 1 shows the error statistics in thecase of the baseline snow depth algorithm (Chang’s algorithm) and the current algorithm. The
current algorithm estimates of depth are compared with measurements from the World
Meteorological Organizations Global Summary of the Day data set (WMO Resolution 40) obtained
from NOAA National Climate Data Center. Seasonal (November through May) root mean squared
error (RMSE) statistics and mean difference between (bias) AMSR-E and in situ measurements(bias) are computed for all available stations in the Northern Hemisphere and then averaged. Results
are stratified for different forest cover fractions and show that with the exception of ff > 0.5, thecurrent AMSR-E depth algorithm has less bias and a lower RMSE compared to the heritage
algorithm, demonstrating that there is an improvement for the new approach.
In a recent study by Tedesco and Narvekar (2010) the NASA algorithm performed comparably
well compared with other snow depth data sets. Furthermore, the standard JAXA snow depth
algorithm was extensively compared with the research algorithms of Koike and another by Palloscia,
both of which were updated in 2010 (JAXA, 2010).

The intercomparison, conducted with 1007

WMO snow depth sites, stratified by land cover classes from ESA’s GlobCover product, and seven

Figure 3 Locations of WMO stations used for testing of snow depth algorithm (Kelly, 2009)
stations located in Yakutsk, Siberia, showed that the standard algorithm’s performance was better
7-13

NDX-120015A

than the other two products for the WMO data and for mapping snow cover extent while the othertwo products did better for the Siberia site. The study showed that the current snow depth algorithm
met the release and standard algorithm target accuracy of less than ±20 cm in 8 of the 15 land covertypes tested and was between ±20 and ±25 cm in 5 categories and greater than 25 cm in 4 categories.
The other two algorithms had only 2 categories with retrieval accuracies better than ±25 cm. The
current algorithm also had the lowest biases in 8 landcover classes with the Koike algorithm having
the lowest in 6 and Palloscia in 1 class. With the initial improvement of the current algorithm overthe heritage Chang et al. (1987) algorithm (Table 1) and the promising performance from the recent
intercomparison study, the standard algorithm has achieved the standard target accuracy of ±20 cm.
Further work is ongoing, however, to achieve the accuracy goal of ±10 cm.

All Data
ff = 0.0
0.0 < ff ≤ 0.5
ff > 0.5

2002 (n=242)
RMSE (cm)
22.1
24.0
21.4
24.6
23.1
24.6
16.9
16.8

2002 (n=242)
Bias (cm)
0.72
6.49
0.78
11.32
2.07
4.38
2.80
0.15

2003 (n=254)
RMSE (cm)
22.5
24.4
24.0
26.2
21.9
24.0
19.2
18.0

2003 (n=254)
Bias (cm)
2.07
6.07
4.57
5.68
0.92
7.05
0.49
1.10

Table 1. Bias and RMSE for the standard AMSR-E algorithm (bold numbers) and the historicalalgorithm Chang et al. (1987) snow depth estimates for 2002-2003 winter (2002) and 2003-2004
winter (2003) for WMO stations in the northern hemisphere. The statistics are reported for: all
stations combined (All Data); stations with no forest cover (ff = 0.0); fractional forest cover greaterthan zero and less than 0.5 (0.0 < ff ≤ 0.5); and forest fractional coverage greater than 0.5 (ff > 0.5). ndenotes the number of WMO stations used in the assessment. (From Kelly, 2009).

References
Armstrong, R.L., A.T.C. Chang, A. Rango and E.G. Josberger (1993) Snow depth and grain size
relationships with relevance for passive microwave studies. Annals of Glaciology, 17: 171-176.
Aschbacher J. (1989) Land surface studies and atmospheric effects by satellite microwave radiometry.
PhD thesis dissertation, University of Innsbruck, Austria.
Biancamaria, S., Mognard, N., Boone, a, Grippa, M., & Josberger, E. (2008). A satellite snow depth
multi-year average derived from SSM/I for the high latitude regions. Remote Sensing ofEnvironment, 112, 2557 - 2568. doi:10.1016/j.rse.2007.12.002
Brown, R., Derksen, C., & Wang, L. (2010). A multi-data set analysis of variability and change in Arctic
spring snow cover extent, 1967–2008. Journal of Geophysical Research, 115(D16).
doi:10.1029/2010JD013975

7-14

NDX-120015A

Chang A.T.C., Foster J.L. & Hall D. K. (1987). Nimbus 7 SMMR derived global snow coverparameters. Annals of Glaciology, 9, 39-44.
Chang, A. T. C., Foster, J. L., & Hall, D. K. (1996). Effects of Forest on the Snow Parameters Derived
From Microwave Measurements During the Boreas Winter Field Campaign. Hydrological
Processes, 10(12), 1565-1574.
Clifford, D. (2010) Global estimates of snow water equivalent from passive microwave instruments:
history, challenges and future developments, International Journal of Remote Sensing, 31: 3707 3726
Derksen, C. (2008). The contribution of AMSR-E 18.7 and 10.7 GHz measurements to improved boreal
forest snow water equivalent retrievals. Remote Sensing of Environment, 112(5), 2701-2710.
doi:10.1016/j.rse.2008.01.001
Derksen, C., A. Walker, & B. Goodison (2005) Evaluation of passive microwave snow water equivalent
retrievals across the boreal forest/tundra transition of western Canada. Remote Sensing ofEnvironment, 96: 315−327.
Derksen, C., P. Toose, A. Rees, L. Wang, M. English, A. Walker and M. Sturm (2010) Development of atundra-specific snow water equivalent retrieval algorithm for satellite passive microwave data.
Remote Sensing of Environment, 114: 1699-1709.
Dewey, K.F. and R.R. Heim (1982) A digital archive of Northern Hemisphere snow cover, November
1966 through 1980, Bulletin of the American Meteorological Society, 11:167-71.
Foster, J.L., A.T.C. Chang and D.K. Hall (1997) Comparison of snow mass estimates from a prototype
passive microwave snow algorithm, a revised algorithm and snow depth climatology. Remote
Sensing of Environment, 62: 132-142.
Gong, G., Cohen, J., Entekhabi, D., & Ge, Y. (2007). Hemispheric-scale climate response to NorthernEurasia land surface characteristics and snow anomalies. Global and Planetary Change, 56(3-4),
359-370. doi:10.1016/j.gloplacha.2006.07.025
Grippa, M., N. Mognard, T. Le Toan, and E.G. Josberger (2004) Siberia snow depth climatology derived
from SSM/I data using a combined dynamic and static algorithm. Remote Sensing of Environment,
93: 30-41.

7-15

NDX-120015A

Grody, N. C., & Basist, a. N. (1996). Global identification of snowcover using SSM/I measurements.
IEEE Transactions on Geoscience and Remote Sensing, 34(1), 237-249.
doi:10.1109/36.481908.
Gunn, G., C. Duguay, C. Derksen, J. Lemmetyinen and P. Toose (2011) Evaluation of the HUT
modified snow emission model over lake ice using airborne passive microwave measurements,
Remote Sensing of Environment, 115: 233-244.
Hansen, M.C., R.S. DeFries, J.R.G. Townshend, M. Carroll, C. Dimiceli, and R. Sohlberg, (2003)
Global Percent Tree Cover at a spatial resolution of 500 meters: first results of the MODIS
continuous fields algorithm. Earth Interactions, 7(10), 15pp.
JAXA (2010) The 4th inter-comparison study of snow depth algorithms for the GCOM/AMSR2,
JAXA/EORC Joint PI Workshop of Global Environment Observation Mission, 6-9 December 2010.
Jordan, R. (1991) A one-dimensional temperature model for a snow cover. CRREL Special Report 91-16,
64pp.
Kelly, R. (2009). The AMSR-E snow depth algorithm: description and initial results. Journal of theRemote Sensing Society of Japan, 29(1), 307-317.
Kelly, R.E.J. and A.T.C. Chang (2003) Development of a passive microwave global snow depth retrieval
algorithm for SSM/I and AMSR-E data. Radio Science, 38: doi : 10.1029/2002RS002648.
Markus, T., Powell, D. C., & Wang, J. R. (2006). Sensitivity of passive microwave snow depth retrievals
to weather effects and snow evolution. IEEE Transactions on Geoscience and Remote Sensing,44(1), 68-77. doi:10.1109/TGRS.2005.860208.
Matzler, C. (1987) Applications of the interaction of microwaves with the natural snow cover. Remote
Sensing Review, 2:259-392.
Metsamaki, S., S. Antilla, H. Markus, and J. Vepsalainen (2005) A feasible method for fractional snow
cover mapping in boreal zone based on a reflectance model. Remote Sensing of Environment, 95:
77-95.
Robinson, D., & Frei, A. (2000). Seasonal Variability of Northern Hemisphere Snow Extent Using
Visible Satellite Data. The Professional Geographer, 52(2), 307-315.
doi:10.1111/0033-0124.00226

7-16

NDX-120015A

Rosenkranz, P.W. (1998) Water vapor microwave continuum absorption: a comparison of
measurements and models. Radio Science, 33: 919-928.
Savoie, M., R.L. Armstrong, M.J. Brodzik and J. Wang (2009) Atmospheric corrections for
improved satellite passive microwave snow cover retrievals over the Tibet Plateau. Remote
Sensing of Environment, 113:2661-2669.
Sun, C.Y, C.M.U. Neale and J.J. McDonnell (1996) Snow wetness estimates of vegetated terrain fromsatellite passive microwave data. Hydrological Processes, 10: 1619-1628
Tedesco M., J. Pulliainen, P. Pampaloni and M. Hallikainen (2004) Artificial neural network based
techniques for the retrieval of swe and snow depth from SSM/I data. Remote Sensing ofEnvironment, 90: 76-85.

Tsang, L., C. T. Chen, A. T. C. Chang, J, Guo and K. H. Ding (2000) Dense Media Radiative
Transfer Theory Based on Quasicrystalline Approximation with Application to Passive
Microwave Remote Sensing of Snow, Radio Sci., 35(3): 731-49.
Ulaby, F.T. and W.H. Stiles (1980) The active and passive microwave response to snow parameters
2: water equivalent of dry snow. Journal of Geophysical Research - Oceans and Atmospheres,85: 1045-1049.
Walker, A.E. and B.E. Goodison (1993) Discrimination of wet snow cover using passive microwave
satellite data. Annals of Glaciology, 17: 307-311.
Wang, J., & Tedesco, M. (2007). Identification of atmospheric influences on the estimation of snow waterequivalent from AMSR-E measurements. Remote Sensing of Environment, 111(2-3), 398-408.
doi:10.1016/j.rse.2006.10.024
Wulder, M. A., Nelson, T. A., Derksen, C., & Seemann, D. (2007). Snow cover variability across central
Canada (1978–2002) derived from satellite passive microwave data. Climatic Change, 82(1-2),
113-130. doi:10.1007/s10584-006-9148-9

7-17

NDX-120015A

Chapter 8.
Description of the GCOM-W1 AMSR2
Soil Moisture Algorithm

Toshio Koike
School of Engineering
The University of Tokyo

8-1

NDX-120015A

Table of Content
1.

Introduction ........................................................................................................................... 2

2.

Algorithm Overview ............................................................................................................. 4

3.

Theoretical Description ......................................................................................................... 4
3.1.

4.

Radiative Transfer Model of the Soil Surface–Vegetation Layer...................................... 4
Algorithm Implementation.................................................................................................... 6

4.1.

Implementation ................................................................................................................. 6

4.1.1.

Indices ....................................................................................................................... 6

4.1.2.

Lookup Tables ........................................................................................................... 7

4.2.

Input/Output parameters ................................................................................................... 7

4.3.

Ancillary data .................................................................................................................... 7

4.3.1.

1.

Vegetation fractional area.......................................................................................... 7

4.4.

Processing flow ................................................................................................................. 8

4.5.

Example output ................................................................................................................. 9

4.6.

Limitations ........................................................................................................................ 95.

Validation Concept .............................................................................................................. 10

6.

References ........................................................................................................................... 12

Introduction
Land surface hydrological quantities have a significant impact on seasonal changes of andinterannual variations in climate through their interactions with the atmosphere. In particular,
variation of the soil moisture content affects the heat balance of the land surface. Information
regarding soil moisture conditions over large regions is important for the understanding, modeling,
and forecasting of climate change. Over short timescales, soil moisture has a strong effect on land
surface evapotranspiration and plays an important role in the redistribution of precipitation thatreaches the land surface into runoff and evapotranspiration. Therefore, frequent observations of soil
moisture are required for successful disaster prevention and water resources management.
Furthermore, soil moisture acts as a water source for vegetation, and monitoring of soil moisture isimportant for understanding ecological processes.
Microwave remote sensing using satellites is an effective method for collecting global
information on land surface hydrological quantities. This method offers two particular advantages: itenables us to periodically perform observations over large regions, regardless of whether it is nightor day, without the microwaves being affected by the atmosphere, and these instruments are highly

8-2

NDX-120015A

sensitive to land surface hydrological parameters, since liquid water has an extremely high dielectric
constant in the microwave band compared with soil.
Investigation of algorithms describing soil moisture began with experimental research using
aircraft-based1) and ground-based2) microwave radiometers. Because vegetation also contains
moisture, soil moisture can only be estimated accurately from satellite observation data by correctly
accounting for the effects of vegetation. Radiative transfer is complicated owing to the wide range ofvegetation properties (particularly with regard to shape, size, and spatial distribution) found in thelong wavelengths from the 1-GHz band to the 6-GHz band; however, radiation and absorption arethought to be the dominant processes, and a good correlation between optical thickness 𝜏𝑐 andvegetation water content 𝑊𝑐 has been demonstrated both experimentally3) and theoretically4). The

relationship between 𝜏𝑐 and 𝑊𝑐 can be expressed as a linear equation, and is used in many soil

moisture algorithms5),6),7) under the assumption of uniform vegetation cover. The microwave

radiometer footprint for satellite observations is several tens of kilometers; however, few land
surface regions globally are uniformly covered with vegetation at this scale. Fujii et al. 8) introduced
vegetation fractional area into their soil moisture algorithm to consider the effects of vegetation
cover; this algorithm provides the basis for the AMSR2 soil moisture algorithm.

8-3

NDX-120015A

2.

Algorithm Overview
Fuiji et al employed the 10–36 GHz algorithm to the AMSR2 L1R product to estimate soil

moisture8) In this algorithm, a lookup table (LUT) method is used for the estimation of soil moisturefrom the L1R brightness temperatures. Because the water content of vegetation affects the sensitivityof the microwave remote sensing of soil moisture, we used a method that simultaneously retrieves
the soil moisture and vegetation water content from two indices, 𝑃𝑃 and 𝐼𝐼𝐼 (the polarization andfrequency differences, respectively, divided by the average value of brightness temperature). Here,the vegetation cover correction of the LUT method was performed using the normalized difference
vegetation index (NDVI) published as part of the moderate resolution imaging spectroradiometer(MODIS) vegetation indices (16-Day L3 Global 1km V5) produced by the land processes distributed
active archive center (LP-DAAC).

3.

Theoretical Description

3.1.

Radiative Transfer Model of the Soil Surface–Vegetation Layer

The low-frequency microwave radiative transfer used in soil moisture calculations is affected
only slightly by the atmosphere. If these effects are ignored, the microwave brightness temperature
𝑇𝑇𝑐,𝑝 of a land surface covered uniformly by vegetation, as observed by satellites, can be expressed

as follows.

𝑇𝑇𝑐,𝑝 = 𝑇𝑇𝑠,𝑝 ∙ 𝑒 −𝜏𝑐 + (1 − 𝜔𝑐 )(1 − 𝑒 −𝜏𝑐 )𝑇𝑐+(1 − 𝜔𝑐 )(1 − 𝑒 −𝜏𝑐 ) ∙ 𝑇𝑐 ∙ 𝛤(𝜃, 𝑝) ∙ 𝑒 −𝜏𝑐(1)where the subscript 𝑝 represents the polarization of the waves; the subscript 𝑐 refers to acomponent related to the vegetation; 𝑇𝑇𝑠,𝑝 is the microwave brightness temperature of the soil

surface; 𝜏𝑐 and 𝜔𝑐 are the optical thickness and single scattering albedo of the vegetation layer,

respectively; 𝑇𝑐 is the physical temperature of the vegetation; 𝛤 is the Fresnel power reflectivity;and 𝜃 is the incidence angle, which is fixed at 55 ° for the AMSR-E. The first term on theright-hand side of Eq. (1) represents the radiation from the ground surface dissipated by thevegetation layer; the second and third terms represent the upward radiation from the vegetation layer
itself and the reflection of the downward radiation, respectively.

In the case of sufficiently moist soil, the microwave brightness temperature 𝑇𝑇𝑠,𝑝 can beapproximated by the product of the soil surface emissivity and the soil physical temperature. Under

thermal equilibrium conditions, this emissivity can be derived from reflectivity. For a smooth surface,
the reflectivity is equal to the Fresnel power reflectivity obtained from the incidence angle and thesoil mixture dielectric constant, where the soil mixture dielectric constant is found by substituting
8-4

NDX-120015A

the soil moisture into Dobson’s model9). The emissivity is also affected by the roughness of the soil
surface, and 𝑇𝑇𝑠,𝑝 is given by the following equation10).

′2𝑇𝑇𝑠,𝑝 = �1 − {(1 − 𝑄) ∙ 𝛤(𝜃, 𝑝) + 𝑄 ∙ 𝛤(𝜃, 𝑞)} ∙ 𝑒 −ℎ cos 𝜃 � ∙ 𝑇𝑠

(2)

where 𝑇𝑠 is the soil physical temperature, 𝑄 is the polarization mixing ratio, ℎ′ is the roughness
height, and 𝑞 is the opposite polarization from 𝑝. 𝑄 and ℎ′ are constants that depend on the′2

surface roughness of the soil. Here, we define 𝐻 = 𝑒 −ℎ 𝑐𝑜𝑠 𝜃 . Values of 𝑄 and 𝐻 for each band
(10 GHz and 36 GHz) can be determined by comparing the brightness temperature data observed by

the AMSR-E with the in situ data at the Mongolia validation site (Table 1).
Only soil surface properties are considered in the soil radiative transfer model given by Eq. (2);however, it has been shown that radiative transfer within the soil layer is important under dry
conditions11). To incorporate this aspect of the soil radiative transfer process into the algorithm, wedeveloped a four-stream fast model for soil that combines the soil radiative transfer model with theconventional soil surface model. In the four-stream fast model, the soil layer is treated as a systemcomposed of an absorbing and radiating medium with a dielectric constant derived from the soil
moisture, embedded with a scattering material composed of many orb-shaped soil particles andlacking an upper boundary. This is combined with a soil surface model that represents the top of thesoil layer, which has no thickness and exhibits only the effects of refraction and reflection. When thesoil is dry, the extinction process has a strong effect in the soil layer. Conversely, as the soil moistureincreases, the effect of the soil layer disappears and the effect of radiation on the soil surface
becomes dominant.
In relation to the vegetation, the optical thickness 𝜏𝑐 and single scattering albedo 𝜔𝑐 in theradiative transfer model are strictly determined using characteristics such as the water content, shape,size, orientation, number, density, and distribution of voids of each of the components that make up
the vegetation (e.g., leaves, stems, and branches). Jackson and Schmugge (1991)3) demonstrated that𝜏𝑐 is highly correlated with the overall water content 𝑊𝑐 (kg/m2) of vegetation and proposed thefollowing linear equation based on experimental results obtained using existing ground-based andaircraft observations.
𝜏𝑐 = 𝑇 ∙ 𝑊𝑐

(3)

where 𝑇 is a vegetation parameter and depends on the type of vegetation. Quantitative information
about 𝜔𝑐 is inadequate because, other than minor scattering at lower frequencies, scattering

depends primarily on the geometrical characteristics of each of the components comprising thevegetation. In this algorithm, fixed values obtained from a field experiment with GBMR are used.
The values of 𝑇 and 𝜔𝑐 used in the algorithm are presented in Table 1.
8-5

NDX-120015A

Table 1: Roughness and vegetation parameters10.65 GHz
V
𝐻

Roughness

𝑄

Vegetation

36.5 GHz
H

V

H

0.873

0.680

0.189

0.344

b

1.93

1.25

3.1

2.9

𝜔𝑐

0.061

0.063

0.061

0.060

The radiative transfer model described above is valid for cases in which the footprint is covered
uniformly with vegetation. According to actual satellite measurements, such uniform vegetation
cover exists in only a few regions globally. In cases where the model domain contains a mixture ofbare and vegetated areas, the microwave brightness temperature 𝑇𝑇𝑝 measured by the satellite canbe expressed by a linear equation obtained from Eqs. (1) and (2):𝑇𝑇𝑝 = (1 − 𝑓𝑐 ) ∙ 𝑇𝑇𝑠,𝑝 + 𝑓𝑐 ∙ 𝑇𝑇𝑐,𝑝(4)

where 𝑓𝑐 is the fractional vegetation cover.

4.

Algorithm Implementation

4.1.

Implementation

4.1.1.

Indices

A radiative transfer model can be used to determine the radiative brightness temperature when allparameters are known, and an inverse analysis method is required for estimating soil moisture using
the observed brightness temperatures. Here, we used a method that allowed simultaneous retrieval of
the soil moisture and vegetation water content from PI2) and ISW12), which represent the polarizationand frequency differences, respectively, divided by the average value of brightness temperature as
defined by Eqs. (5) and (6).
𝑃𝑃 =𝑇𝑇𝑉 − 𝑇𝑇𝐻

1
(𝑇𝑇𝑉 + 𝑇𝑇𝐻 )
2

𝑃𝐼𝑊 =

(5)

𝑇𝑇𝐻,𝑖 − 𝑇𝑇𝐻,𝑗1
�𝑇𝑇𝐻,𝑗 + 𝑇𝑇𝐻,𝑗 �
2(6)

where 𝑇𝑇𝑉 and 𝑇𝑇𝐻 are the microwave radiative brightness temperatures of the vertical and
8-6

NDX-120015A

horizontal polarizations, respectively. The subscripts 𝑖 and 𝑗 in Eq. (6) indicate high and low

frequencies, respectively. There are four unknowns in the radiative transfer model described by Eqs.
(1)–(4) and its related assumptions: soil moisture, vegetation water content, soil physical temperature,and vegetation physical temperature. The PI and ISW depend primarily on soil moisture andvegetation water content and not considerably on the physical temperatures, since the effects of thephysical temperatures decrease as a result of dividing by the average brightness temperature.

4.1.2.

Lookup Tables

It is possible to obtain the soil moisture and vegetation water content corresponding to the PI

and ISW calculated using microwave brightness temperatures: the technique used to correlate these
two indices with the soil moisture and vegetation water content is known as the lookup table method.

Lookup tables are used as inverse analysis tables when retrieving the soil moisture and vegetationwater content from the microwave brightness temperatures and are created in advance for each valueof fractional vegetation cover from the results of radiative transfer value simulations using Eq. (4).
The simulation conditions in this algorithm are as follows:

soil moisture: range of 0.000–0.600 m3/m3, step size of 0.001 m3/m3;vegetation water content: range of 0.000–1.800 kg/m2, step size of 0.001 kg/m2;



fractional vegetation cover: range of 1–100%, step size of 1%; andsoil and vegetation physical temperature: 293 K (fixed).

4.2.

Input/Output Parameters

The input parameters of the soil moisture algorithm include the swath granule brightness
temperature and ancillary data. JAXA AMSR2 L1R swath data form the basis of the soil moistureretrievals and are obtained at EORC as part of the processing chain. The ancillary datasets aredescribed in the next section.
The output parameters of soil moisture algorithm include the volumetric soil moisture andvegetation water contents, although only volumetric soil moisture is saved in the AMSR2 L2 swath
data.

4.3.

Ancillary Data

4.3.1.

Vegetation Fractional Area

We used the NDVI data published as part of the MODIS vegetation indices (16-Day L3 Global 1
km V5) by the LP-DAAC to create a global fractional vegetation cover dataset. MODIS data are

8-7

NDX-120015A

available from two satellites, Terra and Aqua, and the 16-day composite 1-km resolution MODIS
vegetation indices are typically provided as MOD13A2 for Terra and MYD13A2 for Aqua. The
quality flag was checked before data processing to ensure quality and NDVI data excluding cloudy
regions, snow and ice, and water surfaces were used. The fractional vegetation cover was calculated
from NDVI data using the model developed by Carlson and Ripley (1997)13). They performed
simulations using a radiative transfer model for the vegetation–soil–atmosphere, and showed that the
NDVI changes up to total surface cover by vegetation depend mainly on the fractional vegetation
coverage. They then proposed the following fractional vegetation coverage model.

𝑓𝑐 = �

𝑁𝐷𝑉𝑃 − 𝑁𝐷𝑉𝑃𝑚𝑖𝑛 2�
𝑁𝐷𝑉𝑃𝑚𝑎𝑥 − 𝑁𝐷𝑉𝑃𝑚𝑖𝑛(7)

where 𝑁𝐷𝑉𝑃𝑚𝑎𝑥 and 𝑁𝐷𝑉𝑃𝑚𝑖𝑛 are the 𝑁𝐷𝑉𝑃 when the fractional vegetation coverage reaches
100% and the NDVI for bare ground, respectively. In this algorithm, 𝑁𝐷𝑉𝑃𝑚𝑎𝑥 was taken to be 0.7

and 𝑁𝐷𝑉𝑃𝑚𝑖𝑛 was taken to be 0.2 by referring to the variation in 𝑁𝐷𝑉𝑃 for wheat fields in thegrain belt of the state of Victoria in Australia.

Here, two types of processing were carried out to create in advance the fractional vegetation
coverage dataset from the using Eq.(7). The first processing was the spatial averaging. To match the
AMSR-E footprint size, spatial averaging was performed using a 3333 pixel circular boxcar filter at
0.05° grid points. The area of this filter is approximately 797 km2. The AMSR-E footprint is
elliptical and its size is different in each frequency band. However, these factors are not considered
for the dataset in the interests of computational cost in this algorithm. The second processing was
temporal averaging. The MOD13A2 and MYD13A2 both have compositing periods of 16 days, but
the periods start on different days. One starts in the middle of the compositing period of the other.
The temporal averaging was performed every 8 days using the two products to reduce temporaldiscontinuities. The global fractional vegetation coverage dataset was created for the period from
June 2002, when Aqua measurements began, to July 2012.
In the JAXA standard algorithm, nominal values of fractional vegetation coverage, which areaverage for 9 years from June 2002, are used due to limitations of processing system.

4.4.

Processing Flow

Fig. 1 shows an overview of the data processing in this algorithm. First, the fractional vegetation
coverage, 𝑓𝑐 , of each footprint is retrieved using the 𝑓𝑐 dataset and information for the observation

date, latitude and longitude. The appropriate lookup table is then selected. Second, the PI and ISW
are calculated from the 10 and 36 GHz brightness temperature data checked the quality. Finally, thesoil moisture and vegetation water content are retrieved using the lookup table. The ancillary data for

8-8

NDX-120015A

this algorithm only consist of the fractional vegetation coverage information calculated from the
MODIS data and do not require further information.

SATELLITE OBSERVATION
Brightness Temperature (TB)
10GHz (V, H) , 36GHz(V)

Latitude, Longitude,
and Observation date
MODIS NDVI data

Index of Soil Wetness(Koike,1996)
ISW =

1

TB 36 H − TB10 H(TB36 H + TB10 H )2

fractional vegetation coverage
fc (lat, lon, date)

Data source:
Terra/Aqua MODIS 16days, 1km ISIN
Sampling point: 0.05x0.05 grid
Resolution: approx. 33x33 km
fc-model: Carlson and replay, 1997

Polarization Index(Shimonetta,1998)
TB10V − TB10 H(TB10V + TB10 H )Vfc

Look-up Table (LUT)

LUT ( Vfc)

PI

=

1

2

( ISW, PI )
 ( SM, VWC )TBbareSoil Moisture (SM)
Vegetation Water Cont.(VWC)

DATA PROCESSING

Radiative Transfer Model for
vegetated land surface.
TB= (1-fc)TBbare + fcTBvege

bare soil
ANCILLARY DATA

TBvegevegetation

OFF-LINE PROCESSING

Fig. 1: Flowchart of the AMSR2 soil moisture algorithm.

4.5.

Example output

Fig. 2 shows two example outputs extracted from browse data at “the JASMES for water cycle”(http://kuroshio.eorc.jaxa.jp/JASMES/WC.html), which provides information of satellite data related
to the water cycle.
(a) Ascending(b) Descending

Fig. 2: Global map of AMSR2 Soil mositure products in May 20 2013 for (a) acending and (b)descending. The units is volmetric soil mosture content in %.

4.6.
Limitations
There is a tendency to overestimate soil moisture under the extreme dry condition such as

8-9

NDX-120015A

the desert.


Soil moisture around coastline and large lake often shows erroneously high value due towater surface effects.

In the area covered with widely distributed dense vegetation such as Amazon forest andAfrotropical forest, soil moisture shows low value because the signal from soil surface will
be absorbed in the vegetation layer.



Soil moisture also shows low value in heavy precipitation area and mountain area.



RFI in the X-band may be affecting on the SMC product at many locations in the Europe
and Japan.

5.

Validation Concept

There are two concepts for validation of soil moisture products. One is validation at foot print
scale. Other is global applicability.
To get validation data for the first concept, we need to establish field observation network,
considering the radiometer foot print size and land surface heterogeneity. Long-term and high quality
data has been obtained at several sites as shown in the Fig. 3. Mongolia and Khon Kaen (Thailand)sites for soil moisture observation were established in 2000 and 2010. These sites have been
operated by the AMSR/AMSR-E project and were continued today through the AMSR2 projects incooperation with Hirosima University. Yanco (Australia) is also a validation site established in 2012
for soil moisture observation in cooperation with Monash Univeristy. In addition, some validationsites are operated by USGS. Fig. 4 shows a result of comparison of AMSR2 soil moisture within-situ data at Yanco site in Australia from August 2012 to April 2013. In this case, mean absolute
error is about 3.1 % (0.031cm3/cm3) in volumetric soil moisture.
To confirm global applicability, published datasets such as CEOP (http://www.ceop.net) and ISMN
(http://ismn.geo.tuwien.ac.at/) will be useful. However in-situ measurement is point scale in general.
Scale gap between in-situ measurement and AMSR2 footprint may lead to significant uncertainties.
We need to check the land surface heterogeneity around the site carefully to use the point
measurement data.

8-10

NDX-120015A

Mongolia

Mongolia
Khon Kean (Thailand)

Yanco (Australia)

Khon Kean (Thailand)

Yanco (Australia)

JAXA Tower
WeatherstationFig. 3: AMSR2 validation sites

(a)

(b)

2012

2013
*MAE: Mean absolute error

Fig. 4: Comparison of AMSR2 soil moisture with in-situ data at Yanco(Australia). (a) timeseries graph of soil moisture. Red and blue plots are AMSR2 soil moisture in therespective orbits of ascending and descending. Gray plots is in-situ data. Black line isaverage of in-situ data. (b) 2D histogram of AMSR2 soil moisture (Y-axis) versus
in-situ(X-axis).

8-11

NDX-120015A

6.
1)

References
F.T. Ulaby, M.M. Razani, and M.C. Dobson: Effects of Vegetation Cover on the Microwave
Radiometric Sensitivity to Soil Moisture, IEEE Transactions on Geoscience and Remote
Sensing, GE-21(1) pp.51-61, 1983.

2)

S. Paloscia and P. Pampaloni: Microwave Polarization Index for Monitoring Vegetation
Growth, IEEE Transactions on Geoscience and Remote Sensing, 26(5), pp.617-621, 1988.

3)

T. J. Jackson and T. Schmugge: Vegetation Effects on the Microwave Emission of Soils,
Remote Sensing of the Environment,, 36, pp.203-212, 1991.

4)

D. M. L. Vine and M. A. Karam: Dependence of Attenuation in a Vegetation Canopy on
Frequency and Plant Water Content, IEEE Transactions on Geoscience and Remote Sensing,34(5), pp.1090-1096, 1996.

5)

T. J. Jackson: Measuring surface soil moisture using passive microwave remote sensing,Hydrological Processes, 7, pp.139-152. 1993.

6)

E. G. Njoku, T. L. Jackson, V. Lakshmi, T. Chan, and S. V. Nghiem: Soil Moisture Retrieval
from AMSR-E, IEEE Transactions on Geoscience and Remote Sensing, 41(2), pp.215-229,
2003.

7)

T. Koike, Y. Nakamura, I. Kaihotsu, G. Davva, N. Matsuura, K. Tamagawa, and H. Fujii:
Development of an Advanced Microwave Scanning Radiometer (AMSR-E) Algorithm of Soil
Moisture and Vegetation Water Content, Annual Journal of Hydraulic Engineering, JSCE,
48(2), pp.217-222, 2004.

8)

Fujii,H., Koike,T., Imaoka,K.: Improvement of the AMSR-E Algorithm for Soil Moisture
Estimation by Introducing a Fractional Vegetation Coverage Dataset Derived from MODIS
Data, Journal of the remote sensing society of Japan, Vol.29, No.1, pp.282-292., 2009.

9)

M. C. Dobson, F. T. Ulaby, M. T. Hallikainen, and M. A. El-Rayes: Microwave Dielectric
Behavior of Wet Soil - Part II: Dielectric Mixing Models, IEEE Transactions on Geoscienceand Remote Sensing, GE-23(1), pp.35-46, 1985.

10) J. R. Wang and B. J. Choudhury: Remote Sensing of Soil Moisture Content over Bare Field at
1.4 GHz Frequency, Journal of Geophysical Research, 86, pp.5277-5282, 1981.
11) H. Fujii, T. Koike, T. Ohta, H. Ishidaira, T. Jackson, and G. Heathman: Soil moistureobservation under different vegetation conditions by GBMR, Proc. of 2000 IEEE IGARSS,
pp.1268-1270, 2000.

8-12

NDX-120015A

12) T. Koike, T. Tsukamoto, T. Kumakura, M. Lu: Spatial and Seasonal Distribution of Surface
Wetness Derived from Satellite Data, Proceeding of the International Workshop on
Macro-Scale Hydrological Modeling, pp.87-90, 1996.
13) T. N. Carlson and D. A. Ripley: On the relation between NDVI and fractional vegetation
cover and leaf area index, Remote Sensing of the Environment, 62, pp.241-252, 1997.

8-13