# Atmospheric Correction at AERONET Locations: A New Science and Validation Data Set

Yujie Wang^1,2^, Alexei Lyapustin^1,2^, Jeffery L. Privette^3^, Jeffery
T. Morisette^2^, Brent Holben^2^

^1^GEST center, University of Maryland Baltimore county, Catonsville, MD
21228

^2^NASA Goddard Space Flight Center, Greenbelt, MD 20771

^3^NCDC, NOAA Satellite and Information Service, Asheville, NC 28801

**Abstract**

This paper describes an AERONET-based Surface Reflectance Validation
Network (ASRVN) and its dataset of spectral surface bidirectional
reflectance and albedo based on Moderate Resolution Imaging
Spectroradiometer (MODIS) TERRA and AQUA data. The ASRVN is an
operational data collection and processing system. It receives 5050
km^2^ subsets of MODIS L1B data from MODIS adaptive processing system
(MODAPS) and Aerosol Robotic Network (AERONET) aerosol and water vapor
information. Then it performs an atmospheric correction for about 100
AERONET sites based on accurate radiative transfer theory with complex
quality control of the input data. The ASRVN processing software
consists of L1B data gridding algorithm, a new cloud mask algorithm
based on a time series analysis, and an atmospheric correction algorithm
using ancillary AERONET aerosol and water vapor data. The atmospheric
correction is achieved by fitting the MODIS top of atmosphere
measurements, accumulated for 16-day interval, with theoretical
reflectance parameterized in terms of coefficients of the Li-Sparse
Ross-Thick (LSRT) model of the bidirectional reflectance factor (BRF).
The ASRVN takes several steps to ensure high quality of results: 1) a
cloud mask algorithm filters opaque clouds; 2) an aerosol filter has
been developed to filter residual semi-transparent and sub-pixel clouds,
as well as cases with high inhomogeneity of aerosols in the processing
area; 3) imposing requirement of consistency of the new solution with
previously retrieved BRF and albedo; 4) rapid adjustment of the 16-day
retrieval to the surface changes using the last day of measurements; and
5) development of seasonal back-up spectral BRF database to increase
data coverage. The ASRVN provides a gapless or near-gapless coverage for
the processing area. The gaps, caused by clouds, are filled most
naturally with the latest solution for a given pixel. The ASRVN products
include three parameters of LSRT model (*k^L^*, *k^G^*, *k^V^*), surface
albedo, NBRF (a normalized BRF computed for a standard viewing geometry,
VZA=0^o^, SZA=45^o^), and IBRF (instantaneous, or one angle, BRF value
derived from the last day of MODIS measurement for specific viewing
geometry) for MODIS 500m bands 1-7. The results are produced daily at
resolution of 1 km in gridded format. We also provide cloud mask,
quality flag and a browse bitmap image. The ASRVN dataset, including 6
years of MODIS TERRA and 1.5 years of MODIS AQUA data, is available now
as a standard MODIS product (MODASRVN) which can be accessed through
LAADS website (Level 1 and atmosphere archive and Distribution System,
http://ladsweb.nascom.nasa.gov/data/search.html). It can be used for a
wide range of applications including validation analysis and science
research.

## Introduction

Validation of the moderate resolution (\~1km) surface reflectance
products, including spectral bi-directional reflectance factor (BRF),
albedo, vegetation indices and others, is an important component of the
Earth Observing System (EOS) \[1\] and the National Polar Orbiting
Environmental Satellite System (NPOESS) programs. Its goal is to
establish the accuracy of environmental data products on regional and
global scales for a broad range of atmospheric and surface conditions.
The EOS program has developed a multi-level strategy with a strong field
campaign component \[2-5\]. The field measurements required for *direct*
validation analysis provide a detailed and comprehensive look at the
local properties, but they usually involve significant resources, are
subject to weather uncertainties, and are strongly limited in temporal
and spatial coverage. Due to these constraints, recent validation
efforts have proposed that product accuracy assessment should also
utilize a globally representative sample of sites to complement the
direct validation sites \[6\]. This concept has been endorsed by the
Committee on Earth Observing Satellites (CEOS) as the Bench-mark Land
Multisite Analysis and Intercomparision of Products (BELMANIP) \[7\].

In this paper, we present an alternative validation approach for
moderate resolution global surface reflectance products over Aerosol
Robotic Network (AERONET) sunphotometer sites \[8\]. The idea is to
collect the best ancillary information on atmospheric aerosol and water
vapor, and perform an independent atmospheric correction of satellite
measurements based on accurate radiative transfer theory with high
quality control of the input data and results. In the past several
years, we have implemented this idea in the AERONET-based Surface
Reflectance Validation Network (ASRVN) which is an automated data
collection and processing system residing on a dedicated workstation.
ASRVN operationally receives the satellite sensors' L1B data (currently
MODIS TERRA and AQUA from the Goddard's MODAPS and MISR from Langley's
DAAC) and aerosol and water vapor information from the AERONET server.
After a successful test of data integrity and completeness, ASRVN
automatically performs rigorous atmospheric correction on each sensor's
data, creating a sensor-specific record of spectral BRF, albedo and
derivative products (e.g., vegetation index) over more than 100 AERONET
sites globally. For each site, products are stored in a gridded format
at 1 km resolution for an area of 50x50 km^2^ for MODIS and 32x32 km^2^
for MISR. Currently, the ASRVN dataset contains 6 years of MODIS TERRA
data (since 2000), and 1.5 years of MODIS AQUA data. The AQUA dataset
will be completed together with MODIS Collection 5 land re-processing,
which produces the MODIS AQUA subsets for the ASRVN.

Many of the ideas implemented in ASRVN originated in the MODIS
BRDF/Albedo algorithm \[9-10\] similarly based on the time series
analysis within 16-day intervals. In this paper, we will be using the
term BRF (bi-directional reflectance factor) rather than more commonly
used BRDF (bi-directional reflectance distribution function), which is
just a factor of smaller than BRF.

The ASRVN algorithm for MISR data was presented previously \[11\]. This
paper describes processing algorithm for MODIS measurements (sections
2-3). Examples of ASRVN dataset are given in section 4. The paper
concludes with a summary.

## ASRVN Infrastructure

The BRF retrievals from MODIS data use several clear-sky measurements
acquired on successive observations of the same area under different
viewing geometries. The *ASRVN* processing starts by gridding the latest
MODIS swath data (generated operationally over each AERONET site by the
MODAPS production system) and placing the resulting *ASRVN* *Tile* in a
16-day processing Queue, which implements a sliding temporal window
algorithm. Note the *ASRVN Tile* is distinct from the much larger MODIS
Land Team's production Tile. Next, it finds a relevant AERONET aerosol
and water vapor record and computes radiative transfer functions
required for atmospheric correction. After ensuring the quality of input
data by filtering clouds and spatially variable aerosols, the *ASRVN*
algorithm performs atmospheric correction and checks the quality of the
final solution. The next three sections will describe pre-processing
steps preceding the BRF retrieval.

### 2.1 Ancillary Data

AERONET sunphotometers sample the direct solar radiation each 15
minutes, and sample diffuse sky radiance over a wide range of angles
every 60 minutes during the daytime. AERONET's automated processing
system generates AOT and column water vapor from the direct solar
measurements. Typical AOT uncertainty for a field instrument is
0.01-0.02 and is spectrally dependent. The inversion algorithm \[12\]
uses almucantar sky measurements to retrieve aerosol microphysical
properties (particle size distribution and refractive index) and
concentration. AERONET applies several tests to ensure reliability of
retrievals, such as SZA 45^o^, AOT~0.44~ 0.4, and that there were at
least 21 independent angles used in each inversion. The tests analyze
sensitivity of retrievals to the single scattering albedo, and to the
phase function at large scattering angles. These quality assurance tests
significantly reduce the number of complete aerosol characterization
records, as compared to the number of AOT records.

The *ASRVN* algorithm starts with the selection of AERONET aerosol
optical thickness and column water vapor values within 30 minutes of
satellite overpass. If these conditions are met, the algorithm selects
the inversion record with aerosol microphysical parameters and size
distribution within 2 hours of overpass. Otherwise, it uses an aerosol
climatology model for a given location derived from multi-year AERONET
statistics of reliable retrievals \[13\]. Because full AERONET
inversions are less accurate when aerosol concentration is low, the
climatology background aerosol model is always used in our algorithm for
clear atmospheric conditions (currently defined as AOT~0.44~0.3). Our
testing and earlier prototyping with MISR data \[11\] demonstrated that
the aerosol climatology significantly improves stability of the time
series of derived surface albedo.

Following the selection of aerosol parameters, *ASRVN* algorithm
calculates the aerosol optical thickness (AOT~MIE~), single scattering
albedo, and scattering phase function using a look-up table approach
\[14\]. Depending on AERONET sphericity index, either a spherical (Mie)
aerosol model or a model of spheroids is used. The above calculations
provide the spectral dependence of extinction (AOT~MIE~) in the MODIS
wavelengths. However, the AOT from direct solar measurements may differ
from AOT~MIE~ for a number of reasons, from the time difference between
inversion and direct AOT measurement to uncertainties associated with
the AERONET inversion algorithm \[15\] etc. For this reason, AOT~MIE~ is
further scaled by fitting it to the measured AOT at three AERONET
wavelengths (0.44, 0.67, and 0.87 m). Once aerosol optical parameters
are defined, the radiative transfer model SHARM \[16\] calculates the
required radiative transfer functions for the specific water vapor and
spectral response functions of MODIS TERRA or AQUA instrument using the
Interpolation and Profile Correction (IPC) method \[17\]. A pixel-wise
correction for variations of surface elevation (atmospheric pressure) is
performed using a spectral interpolation method \[18\].

### 2.2 Implementation of Time Series Processing

To execute time series processing (sliding window algorithm), *ASRVN*
first grids MODIS level 1B (L1B) calibrated and geolocated data to a
regular 1 km grid. We use the MODIS land gridding algorithm \[19\] with
minor modifications that allow us to better preserve the angular
anisotropy of signals in the gridded data when measured reflectance is
high, for example over snow, thick clouds or water with glint. Next,
gridded MODIS data (*Tiles*) are placed in the processing Queue, which
can hold up to 16 days of successive measurements. The *ASRVN*
processing uses both individual grid cells, also called pixels below,
and fixed-size (2525 km^2^) areas, or blocks, required by the cloud mask
algorithm. In order to organize such processing, we developed a
framework of C++ classes and structures (algorithm-specific Containers).
The class functions are designed to handle processing in the various
time-space scales, for example at the pixel- or block-level, and for a
single (last) day of measurements or all available days in the Queue, or
for a subset of days which satisfy certain requirements (filters). The
data storage in the Queue is efficiently organized using pointers, which
avoids physically moving the previous data in memory when the new data
arrive.

The structure of the Queue is shown schematically in Figure 1. For every
day of observations, MODIS measurements are stored as Layers for
reflective bands 1-7 and thermal bands 31-32, all of which are required
by the cloud mask (CM) algorithm. Besides storing gridded MODIS data
(*Tiles*), the Queue has a dedicated memory (Q-memory) which accumulates
ancillary information about every block and pixel of the surface for the
cloud mask algorithm (*Refcm* data structure). It also keeps information
related to the history of previous retrievals, such as surface BRF
parameters and albedo. Given the daily rate of MODIS observations, the
land surface is relatively static over most 16-day periods. Therefore,
knowledge of the previous surface state significantly enhances both the
accuracy of the cloud detection and the quality of atmospheric
correction by imposing a requirement of consistency of the time series
of BRF and albedo.

### 2.3 Data Quality Control: Cloud Mask and Aerosol Filter

From the start, *ASRVN* was designed to work with a multipixel area
rather than a single pixel centered at an AERONET location in order to
provide optimal visual control over the input and output results. Visual
analysis of RGB images is superior for complex data quality assessment
and troubleshooting situations which is rarely achievable with a pixel-
(or point-) level analysis.

Although AERONET produces an internal cloud mask \[20\], it is not
sufficient for atmospheric correction over the 5050 km^2^ processing
area. For example, the sunphotometer may provide AOT measurements from a
direct sun view through a gap in the clouds. Usually, the
sunphotometer's time of measurement differs from satellite overpass
time, in which case changes in cloudiness may have occurred. For these
reasons, we implemented a new *MAIAC* Cloud Mask algorithm \[21\] as
part of the *ASRVN* processing. This algorithm uses the time series
analysis and an image-based rather than pixel-based analysis. It
explicitly addresses the problem of cloud searching by identifying a
clear-skies comparison target. *MAIAC* CM constructs the reference
clear-skies image of the surface and stores it along with other
ancillary information about reflectance and brightness temperature for
every surface block (2525 km^2^). This ancillary information, required
for cloud detection, is continuously updated with the latest cloud-free
measurements. This allows the reference image to dynamically adjust to
gradual (seasonal) and rapid surface changes caused by snowfall, fires,
floods, etc. The CM algorithm is enhanced by an internal dynamic
land-water-snow classification, which allows processing flexibility over
varying surface types.

Our experience with MODIS data processing dictated the need for
additional data screening. This is required when aerosols have
significant spatial variation and a single AERONET AOT value does not
represent the full processing area, as well as in cases of undetected,
usually semi-transparent and sub-pixel, clouds. This screening was
implemented through an "aerosol filter". Using the known surface BRF
from previous retrievals, the algorithm computes the pixel-level AOT in
the blue band from the latest MODIS measurements. The AOT retrieval is a
fast algorithm based on a look-up table pre-calculated for a standard
continental aerosol model. The computed AOT is used solely to assess
spatial homogeneity of aerosols over the processing area, and to find
deviations, which usually indicate previously undetected clouds and
sometimes spatially varying aerosols. Specifically, the algorithm
generates an AOT histogram from the non-cloudy pixels, filters the
highest 20% values as possibly cloudy, and finds the average value
(AOT~av~) for the remaining 80% of the pixels. This AOT~av~ is assumed
to represent the average clear-skies aerosol loading over the processing
area, which should correspond to the AERONET AOT. The AOT~av~ is used
next to further filter "suspicious" (contaminated) data as follows: if
the atmosphere is clear (AOT~av~\<0.25), then the algorithm filters only
pixels with the high AOT values exceeding the average by 0.15 or more.
This threshold was determined through trial and error. Otherwise, it
filters the high and low values symmetrically if the difference with
AOT~av~ exceeds 0.15. This relatively simple technique allows us to
filter sub-pixel clouds, contrails and other forms of thin cirrus and
semi-transparent clouds. With the introduction of this additional
filter, we witnessed a dramatic enhancement in the quality of the
*ASRVN* atmospheric correction.

Figure 2 shows examples of the cloud mask algorithm and aerosol filter
for the Goddard Space Flight Center (Greenbelt, Maryland, USA) site. It
shows that the CM algorithm captures most of the opaque clouds, whereas
the aerosol filter captured additional sub-pixel and semi-transparent
clouds and cases of spatially variable aerosols. The example at the
bottom of Figure 2 also shows the detection of cloud shadows by *MAIAC*
CM algorithm. Shadows are detected with a simple threshold algorithm
which compares the latest MODIS measurement (*^meas^*) with predicted
top-of-atmosphere (TOA) reflectance from the previously retrieved BRF
model parameters (*^pred^*):

IF *^meas^* \< *^pred^* - 0.12 *CLOUD_SHADOW*.

Here, we use MODIS wavelength 1.24 m (band 5) which experiences minimal
atmospheric distortions and is usually bright over land so the change of
reflectance due to cloud shadow can be easily detected well above the
sensor noise level.

## Atmospheric Correction

Once the cloud mask, enhanced by aerosol filter, is applied, the *ASRVN*
algorithm filters the time series of MODIS measurements for every single
pixel, and places the remaining good data in a "container". The
container stores measurements along with computed RT functions for the
cloud-free days of the Queue. If the number of good measurements exceeds
3 for a given pixel (see section 3.2), then the coefficients of
Li-Sparse Ross-Thick (LSRT) BRF model \[22\] are computed. The LSRT
model is used in the MODIS BRDF/albedo algorithm \[9-10\].

**3.1 Inversion for LSRT Coefficients**

In the current operational MODIS land processing, the BRF is determined
in two steps: first, the atmospheric correction algorithm derives
surface reflectance for a given observation geometry using a Lambertian
approximation \[23\]. Next, three LSRT coefficients are retrieved from
the time series of surface reflectance accumulated for a 16-day period
\[9\]. The Lambertian assumption simplifies the atmospheric correction
but imparts biases which depend on observation geometry and atmospheric
opacity. Tests show that the Lambertian assumption leads to a more
Lambertian BRF shape while the true BRF is more anisotropic \[24\].

The *ASRVN* algorithm derives surface LSRT coefficients directly by
fitting the radiative transfer solution to the measured TOA reflectance
accumulated over a 4-16 day period. The inversion is based on a high
accuracy semi-analytical Green's function solution \[25-26\], which in
combination with LSRT BRF model provides an explicit parameterization of
TOA reflectance in terms of the surface BRF model parameters
$K$={*k^L^*, *k^G^*, *k^V^*}*^T^*. According to the derivation provided
in the Appendix, the TOA reflectance can be expressed as:

$R(\mu_{0},\mu,\phi) = R^{D}(\mu_{0},\mu,\phi)$+$k^{L}F^{L}(\mu_{0},\mu)$+$k^{G}F^{G}(\mu_{0},\mu,\phi)$+$k^{V}F^{V}(\mu_{0},\mu,\phi)$+$R^{\text{nl}}(\mu_{0},\mu)$,
(1)

where $R^{D}$ is the atmospheric (path) reflectance, and $R^{\text{nl}}$
is a small non-linear term proportional to the product of surface albedo
and spherical albedo of the atmosphere
($R^{\text{nl}} \propto \text{qc}_{0}$). Functions $F^{L}$, $F^{V}$ and
$F^{G}$ depend on geometry and atmospheric conditions. They are weakly
non-linear in *k*-coefficients through the multiple reflection factor
$\alpha = (1 - \text{qc}_{0})^{- 1}\text{.}$

The quasi-linear form of equation (1) leads to a very efficient
iterative minimization algorithm:

*RMSE* =
$\sum_{j}^{}{(r_{j}^{(n)} - F_{j}^{L}k^{L(n)} - F_{j}^{V}k^{V(n)} - F_{j}^{G}k^{G(n)})^{2}}$
=$\underset{\left\{ \overrightarrow{K} \right\}}{\text{min}}$,
$r^{(n)} = R - R^{D} - R^{\text{nl}(n - 1)}$, (2)

where index *j* denotes measurements for different days, and *n* is the
iteration number. Equation (2) provides an explicit least-squares
solution for the kernel weights. In matrix form, the solution is written
as:

${\overrightarrow{K}}^{(n)} = A^{- 1}{\overrightarrow{b}}^{(n)}$, (3)

where

$A = \begin{bmatrix}
\sum_{}^{}{(F_{j}^{L})^{2}} & \sum_{j}^{}{F_{j}^{G}F_{j}^{L}} & \sum_{j}^{}{F_{j}^{V}F_{j}^{L}} \\
\sum_{j}^{}{F_{j}^{G}F_{j}^{L}} & \sum_{j}^{}{(F_{j}^{G})^{2}} & \sum_{j}^{}{F_{j}^{V}F_{j}^{G}} \\
\sum_{j}^{}{F_{j}^{V}F_{j}^{L}} & \sum_{j}^{}{F_{j}^{V}F_{j}^{G}} & \sum_{j}^{}{(F_{j}^{V})^{2}}
\end{bmatrix}$, ${\overrightarrow{b}}^{(n)} = \begin{bmatrix}
\sum_{j}^{}{r_{j}^{(n)}F_{j}^{L}} \\
\sum_{j}^{}{r_{j}^{(n)}F_{j}^{G}} \\
\sum_{j}^{}{r_{j}^{(n)}F_{j}^{V}}
\end{bmatrix}$.

In the first iteration, the small non-linear term is set to zero,
$R_{j}^{\text{nl}(0)} = 0$, and the multiple reflection factor (see
Appendix) is set to one, $\alpha^{(0)} = 1$. These parameters are
updated once after the BRF coefficients are calculated in the first
iteration. Except for snow-covered surfaces, the problem converges with
high accuracy in two iterations because the non-linear terms are small.
Currently, the ASRVN algorithm does not make retrievals over snow.

Prior to inversion, the algorithm checks if the dataset has a sufficient
angular sampling. The MODIS operational BRDF/albedo algorithm \[9\]
makes an inversion if at least 7 cloud-free observations are available
during 16-day period. We studied this problem experimentally using MODIS
data for a number of AERONET sites, varying the minimal required number
of measurements (from 3 to 10) and testing different metrics of angular
sampling. One metric used the magnitude of determinant of the inverse
matrix *A* which shows how different the sampling angles are. Although
such analysis is, perhaps, most straightforward theoretically, we found
it often too restrictive. In the end, a simple criterion was chosen
based on the range of cosine of the view zenith angle (~max~ - ~min~
0.2), which ensures robust and consistent retrievals.

The described algorithm (Eqs. 1-3) has a high computational efficiency.
Compared to the radiative transfer computations, the time required to
evaluate functions $F^{m}$ (*m*=*L*, *V*, *G*) and $R^{\text{nl}}$ is
negligible. The integrals required for these functions (see Appendix)
need to be calculated only once regardless of the number of iterations.
Finally, small variations of the viewing geometry across the processing
area are neglected and the RT calculation is done once per observation
using the geometry of the central pixel of the subset.

### 

### 3.2 Solution Selection and Update

Although the LSRT model leads to an efficient BRF retrieval algorithm,
there are several caveats associated with this model. The LSRT kernels
are not orthogonal, are not positive-only functions, and are normalized
in a somewhat arbitrary fashion that is not linked to radiative transfer
theory. These factors reduce the stability and uniqueness of the
solutions, such that small perturbations in measurements may lead to
significantly different solutions. The high goodness-of-fit at the
measurement angles does not guarantee the correct shape of the retrieved
BRF, and may result in negative BRF values at other angles. The albedo,
being an integral function of BRF, is especially sensitive to an
incorrect BRF shape. For these reasons, we developed several tests to
remove unrealistic solutions.

The initial validation of the solution (see Figure 3) checks that the
maximal difference over all days of the Queue between measured and
computed TOA reflectance does not exceed a specified threshold
($\left| R^{\text{Meas}} - R^{\text{LSRT}} \right|$\>0.08). The day
(measurement) with the highest deviation is excluded from the Queue and
the inversion is repeated. If the number of measurements goes below four
after the exclusion, no retrieval will be made for this pixel.

If a solution provides a good agreement with measurements for all days,
the algorithm verifies that values of the direct-beam albedo (*q*) at
SZA=15^o^, 45^o^, 60^o^ are positive. Finally, the new solution must be
consistent with the previous solution:
$\left| q(\text{45}^{o}) - q^{\text{Pr}\text{ev}}(\text{45}^{o}) \right|$\<(),
where is the band-dependent threshold currently equal to 0.04 (blue),
0.05 (green and red), 0.1 (for spectral region of 0.8 -- 1.6 m) and 0.05
for the shortwave infrared band (2.1 m). Consistency of the time series
of BRF and albedo is characterized by a *status* index. Initially, the
confidence in the solution is low (*status*=0). Each time the new
retrieval agrees with the previous retrieval, *status* increases by 1.
When *status*3, the retrieval is considered reliable.

The thresholds (0.08 and ()) in the LSRT inversion routine are selected,
on one hand, tight enough to reject most of undetected clouds, which
remain the dominant source of errors, and sufficiently loose, on the
other hand, for the solution to adapt to the surface change. The most
pervasive type of change is seasonal variations, related to the spring
green-up and fall senescence at northern latitudes, or greenness
variations caused by wet and dry seasons in tropics. The total seasonal
variation of reflectance (e.g., see Figure 5) is about several absolute
percent in the visible bands, and is significantly higher in the near
infrared (\~0.1-0.2). The ASRVN thresholds for the daily variation are
selected accordingly, and our analysis of a large volume of processed
MODIS data confirms that the ASRVN algorithm does not reject
measurements when surface is changing, even in the agricultural regions
characterized by a rapid reflectance change during harvesting.

When the new solution is validated, the coefficients of the BRF model
and direct-beam albedo *q*(45^o^), stored in the Q-memory, are updated.
The update is done with relaxation, designed to mitigate random noise of
retrievals:

${\overrightarrow{K}}_{\lambda}^{\text{New}}$=
(${\overrightarrow{K}}_{\lambda}^{\text{New}}$ +
${\overrightarrow{K}}_{\lambda}^{\text{Pr}\text{ev}}$)/2. (4)

This method of update increases the quality of the BRF and albedo
product when the surface is relatively stable, but it delays the
response of the solution to surface changes.

Often, the solution for some pixels or the full area cannot be produced
because of lack of clear-skies measurements. In these cases, we assume
that the surface does not change and we fill-in data gaps with the
previous solution for up to a 32-day period. In most cases, this
assumption of a stable surface is reasonable. The gap-filled pixel is
marked as "Extended" in the quality assurance (QA) value, with parameter
QA.*nDelay* giving the number of days since the last update (see section
3.6).

### 3.3 *ASRVN* Products

The *ASRVN* computes two main products at 1 km resolution for seven 500m
MODIS bands, the set of BRF coefficients, and the surface albedo. The
albedo is defined by Equation (A-5a) as a ratio of surface-reflected to
incident radiative fluxes. Thus, it represents a true albedo at a given
solar zenith angle in ambient atmospheric conditions, the value, which
can be directly compared to ground-based measurements.

*ASRVN* also computes several derivative products useful for science
data analysis and validation:

1\) *NBRF* - a BRF Normalized to the common geometry of nadir view and
SZA = 45^o^. This product is analogous to MODIS NBAR (nadir BRF-adjusted
reflectance) product (part of the MOD43 standard product suite). With
the geometry variations removed, the time series of *NBRF* is useful for
studying vegetation phenology, performing surface classification, etc.

2\) *IBRF* -- an instantaneous (or one-angle) BRF value for the specific
viewing geometry of the last day of observations. In essence, *IBRF* is
a reflectance which would be measured if the atmosphere were absent.
This product is calculated from the latest MODIS measurement assuming
that the shape of BRF, known from previous retrievals, didn't change. To
illustrate computation of *IBRF*, we re-write equation (1) for the
measured TOA reflectance as follows:

$R(\mu_{0},\mu,\phi) = R^{D}(\mu_{0},\mu,\phi)$ +
$\text{bR}^{\text{Surf}}(\mu_{0},\mu,\phi)$, (5)

where *R^Surf^* combines all surface related terms and can be calculated
using the previous solution for BRF (*BRF*) and AERONET aerosol
information. *b* is spectrally-dependent scaling factor. Then,

$\text{IBRF}_{\lambda}(\mu_{0},\mu,\phi) = b_{\lambda}\text{BRF}_{\lambda}(\mu_{0},\mu,\phi)$.
(6)

This algorithm (Equation 6) will be referred to as scaling. This
description was given for the purpose of illustration. In reality,
*R^Surf^* is a non-linear function so computing parameter *b* and *IBRF*
is done accurately using the formulas given in the Appendix.

The algorithm computing scaling coefficient (and *IBRF*) is shown in
Figure 3 on the right. First, the algorithm filters measurements which
differ from the theoretically predicted TOA reflectance based on the
previous solution ($R_{Q}^{\text{LSRT}}$) by more than factor of ().
Then, the scaling coefficients are computed, and the consistency
requirement is verified as: 0.5\<*b*\<2 in the usually dark visible
bands, and 0.8\<*b*\<1.2 in the bright near-infrared band. The scaling
coefficient's range of variability is roughly consistent with the
thresholds for the albedo variations discussed in the previous section.
If all conditions are satisfied and the *status* of the pixel is high
(*status*3), then the LSRT BRF parameters of the pixel stored in the
Q-memory are updated with the scaled solution:

${\overrightarrow{K}}_{\lambda}^{\text{New}}$=$\frac{b_{\lambda} + 1}{2}$${\overrightarrow{K}}_{\lambda}^{\text{Pr}\text{ev}}$.
(7)

Based on its definition, *IBRF* is well suited for validation of surface
reflectance product of the operational atmospheric correction algorithm
(standard product MOD09 \[23\]).

The cloud mask, the RGB browse images for MODIS TOA reflectance and for
the ASRVN NBRF and IBRF, and the QA flag, are also standard parts of the
ASRVN product suite. The products are saved in HDF-EOS format files,
which can automatically keep geolocation information and allow data to
be ported to virtually any computer platform, regardless of the byte
order of the native platform. A list of the ASRVN products and their
MODIS products counterparts are given in Table 1, and the definition of
the cloud mask values is provided in Table 2.

### 3.4 Update in Case of Rapid Surface Change

Time series processing is intrinsically controversial when a surface
changes rapidly. Since inversions are generally ill-posed problems, one
desires all available cloud-free measurements and a maximal time window
in order to reduce the *RMSE.* This approach, which reduces the impact
of noise in the data (including that of gridding and of residual clouds)
and which ensures a more robust BRF shape, is best when the surface is
stable throughout the accumulation period. For example, this is the case
for natural ecosystems in mid-latitudinal summers. In contrast,
detecting and tracking surface changes like spring green-up,
agricultural harvesting or fall senescence requires the least possible
number of days in the inversion Queue. Such retrievals may have a
considerable amount of spatial and spectral noise. Indeed, it is
difficult to assess the reliability of solutions when the surface
reflectance is changing rapidly or abruptly, especially given the
possibility of data gaps due to clouds.

While the response of the 16-day solution may be delayed, the *IBRF*
tracks spectral changes immediately. The update of Q-memory with the
latest measurements by equation (7) was found to significantly
accelerate the response of LSRT coefficients
(${\overrightarrow{K}}_{\lambda}$), and hence of the NBRF, to changing
surface conditions. Yet, in some cases we found this insufficient. To
amend results in these cases, we added a separate update of the
Q-solution with *IBRF* based on change detection. This path is shown at
the bottom of Figure 3. Usually, seasonal surface changes related to the
vegetation cycle are accompanied by a correlative change in the red-NIR
bands, for example, a simultaneous decrease of red and increase of NIR
reflectance during green-up, and the opposite changes during fall
senescence and defoliation for broadleaf forests. In our case, the
change detection is based on the top-of-canopy Normalized Difference
Vegetation Index (NDVI) \[27\]. The scene-average NDVI is calculated
using the *IBRF* (NDVI^I^) and then is normalized to a standard viewing
geometry of NBRF. It is compared to the scene-average NDVI^Q^ calculated
using the NBRF stored in the Q-memory (previous reliable solution). A
change is defined as when the difference between the two values exceeds
0.01. At the same time, the red and NIR reflectance should change
accordingly. For example, the following set of rules defines the change
during spring green-up:

NDVI^I^ - NDVI^Q^ \> 0.01,
$\rho_{\text{Re}d}^{I} - \rho_{\text{Re}d}^{Q} < 0$,
$\rho_{\text{NIR}}^{I} - \rho_{\text{NIR}}^{Q} > 0$.

Such an approach filters fluctuations in NDVI caused by factors
unrelated to the vegetation signal, such as residual cloudiness or
undetected cloud shadows. Once the change is detected, the Q-memory is
updated using Eq. (7). Empirically, we found that this approach, which
is responsible for 10-15% of all updates, performs robustly at different
global AERONET locations, for example in North America, Eurasia, or
Africa.

The NBRF and albedo are also updated following an update of
${\overrightarrow{K}}_{\lambda}$ coefficients. This method, which
requires only one last day of measurements when the BRF shape is
reliably known, is similar to the scheme used in the MODIS BRDF/albedo
algorithm \[10\] and its reliability has been tested with a global
composite AVHRR data set \[28\]. An advantage of this strategy is that
once initialized, the algorithm provides near gapless coverage, and a
continuous time series of NBRF and albedo.

Besides providing a faster response to surface change, the developed
update strategy assures fast removal of retrieval artifacts, mainly
residual clouds, which is the most common problem.

### 3.5 Seasonal Ancillary BRF

When no reliable retrievals are made during the past 32 days, which is
usually caused by high cloudiness, the previous retrievals are
considered unreliable and the Q-memory is refreshed with fill values.
After that, the algorithm may take a considerable amount of time to
re-initialize, during which time no results will be produced. To remedy
this situation, we developed the historical 1 km resolution BRF database
from 5 year retrievals, which is used for BRF-scaling when the Q-memory
is being refreshed. This idea of historical back-up BRF, providing the
angular shape of function, was initially proposed and implemented in the
MODIS BRDF/albedo algorithm \[10\]. The ASRVN BRF database contains one
set of spectral K-coefficients for every pixel for each of the four
calendar seasons. Initially, the database is built from a multi-year run
of ASRVN. The images used for averaging, are selected according to the
standard deviation of NBRF (~NBRF~) computed for the processing area.
Empirically, we found that with our data quality assurance, the top 20%
of images with the highest standard deviation usually contain artifacts
from residual clouds or unreliable BRF solutions, whereas the remaining
80% of retrievals have a good quality. Thus, for each site and every
season, we first generate the histogram of ~NBRF~, find the 80%
threshold, and average the images with -values lower than the threshold.

Once the seasonal BRF database is created, it is supported by an offline
background algorithm which updates it with the latest good quality
solution.

### 3.6 Quality Assurance (QA) Flag

For each execution, the algorithm creates a pixel-level QA flag to
indicate the overall quality and the internal processing path. The QA
information consists of 16-bit compound bit-fields, as summarized in
Table 3.

The *QA.overall* field indicates the overall data quality. Four values
are possible:

1)  Good. This means that the *status* of solution is 3 or greater. In
    other words, at least three of the last LSRT retrievals agreed, and
    the calculated IBRF was found to be consistent with the model
    predictions;

2)  Acceptable. In this level, the *status* is not high (\<=2), but the
    calculated IBRF agrees with the model prediction. This case may
    represent a good solution where only a few retrievals are available
    because of a gap in AERONET data or clouds;

3)  Extended. In this case, either the solution was not produced because
    of a lack of clear-sky measurements, or the calculated IBRF did not
    agree well with the model prediction. In this case, the previous
    reliable solution for a given pixel is used to fill-in the values of
    the LSRT coefficients and NBRF. The IBRF in this case is not
    produced;

4)  Not created. This usually happens at the beginning of processing
    when the Q-memory is not yet initialized.

The field *QA.nDelay* gives the number of days since the last update of
the Q-memory. If *QA.nDelay* =0, then this pixel contains the most
recent retrieval. If the Q-memory was not updated for 32 or more days,
the information for a given pixel will be overwritten with the fill
value. If the solution is calculated with (Eq. 7), the *QA.scale* field
will be set to 1 to indicate that the solution is "scaled" from the
previous reliable retrieval. The value of the cloud mask for a given
pixel is stored in the field *QA.cloud*. We also mask pixels which are
adjacent to cloudy pixels, where greater errors in the atmospheric
correction are expected. For such pixels, the *QA.adjCloud* field is set
to 1 (the default value is 0).

## 4. Processing Examples

We selected three AERONET sites, Goddard Space Flight Center (GSFC,
38.9925^o^N, 76.84^o^ W), Mongu, Zambia (15.25^o^S, 23.15^o^E) and Solar
Village, Saudi Arabia (24.91^o^N, 46.41^o^E) to illustrate the ASRVN
dataset. These sites have very different land cover types, atmospheric
conditions and seasonal variations.

-   The GSFC site is located in a northern suburb of Washington, DC. It
    is a mixture of urban residential area, small deciduous broad leaf
    forest stands and small patches of agricultural cropland.

-   The Mongu site is located on the eastern side of Zambezi River in
    western Zambia. The western part of the area is a floodplain mostly
    covered with grasses, and the eastern part is mainly a sandy soil
    with sparse vegetation. Significant biomass burning takes place in
    August and September. In the wet season (November-March), this area
    has a high cloudiness.

-   The Solar Village site is a desert area with stable surface
    conditions. The dominant aerosol source is dust.

Below, we provide several examples to demonstrate *ASRVN* products, data
quality and potential applications. For this study, we used MODAPS
Collection 5 data exclusively.

### 4.1 NBRF time series

The NBRF product has been corrected for both atmospheric effects and
variations of view geometry. Thus, NBRF variations should be closely
related to changes in surface conditions. Figure 4 shows the seasonal
dynamics of the NBRF images for the three sites. The "true color" images
are composed with equal weights from the red, green and blue spectral
bands. The columns show the gridded TOA reflectance, NBRF and the new
cloud mask (the CM legend is given in the figure's caption). Seasonal
changes are easy to observe in the first two columns. The vegetative
cover reaches a maximum in July-August for the GSFC site, whereas at
Mongu, green vegetation is most abundant in January-February, at the
height of the wet season. Surface reflectance at the Solar Village site
exhibits little variation throughout the year as is expected for a
desert location with little vegetation.

To demonstrate the algorithm's performance with different surface types,
we selected two pixels for GSFC that differ in the amount of vegetation
-- a "green" pixel \[pixel (16, 36)\] in the middle of a small deciduous
forest stand to the north-east of the center, and a relatively bright
"urban" pixel \[pixel (46, 3)\] representing a typical residential area
with mixed vegetation at the lower-left corner of the image. For
comparison, we also selected a "bright" pixel \[pixel (10, 20)\] in the
desert region of Solar Village. The locations of these pixels are
indicated by color circles in Fig. 4. The NBRF time series for these
pixels are shown in Fig. 5. For the "green" pixel, the NBRF in the NIR
band increases rapidly during springtime, while the red, green, and blue
reflectances decrease. There is an interesting dynamic between the red
and green signals. In the early spring, reflectance in the red channel
is greater than in the green channel. This is typical of most soils.
With the spring green-up, the red band reflectance decreases
significantly due to chlorophyll absorption while the change in the
green band reflectance is much smaller. During the autumn season, the
bands change in the reverse direction as expected for senescing
vegetation (Fig. 5a). The "urban" pixel shows a similar temporal trend,
but with much smaller amplitude (Fig. 5b). The time series of the
"bright" pixel at Solar Village does not show much variation throughout
the year. The NBRF of band 7 and the NIR, red, green and blue bands
remains around 0.5, 0.42, 0.35, 0.24, and 0.14 respectively. The
variation is about +/- 0.03-0.05 in each band.

The data gap in Fig. 5 in year 2004 is due to the incompleteness of the
MODIS dataset we acquired to date. The remaining gaps will be filled in
upon completion of the MODIS land Collection 5 reprocessing, which also
generates the MODIS subsets for the AERONET sites for the ASRVN.

### 4.2 IBRF vs. NBRF

As discussed in section 3.4, the NBRF, which is retrieved from 4 to 16
days of measurements, may have a delayed response to surface changes. In
contrast, the IBRF, derived from the last day of measurements, tracks
surface spectral changes immediately. This improvement in temporal
sensitivity is sometimes achieved with compromised data quality.

Figure 6 shows a comparison of NBRF and IBRF for the GSFC site for the
spring green-up period (days 95-122) of 2005. There is a 7-day gap in
the AERONET records due to cloudiness after day 101. As a consequence,
the NBRF and IBRF images show a noticeable color difference on day 108,
representing a delayed response of the NBRF to surface change. With the
accumulation of six additional measurements, the color of the NBRF and
IBRF images become consistent again on day 117. The lag in NBRF depends
on the number of available clear-skies measurements and the rate of
surface change.

Note, that the brightness of the IBRF images in Fig. 6 vary with the
view zenith angle. This occurs because the IBRF images are not
normalized, whereas the NBRF images show a stable green-up trend as
expected. This confirms that IBRF is more useful for detecting surface
changes as well as effects of storms, fires, etc. The more stable NBRF
is more suitable for characterizing stable periods, detecting the
long-term trends, and support process models, which require low noise of
the input data.

To illustrate differences and similarities between NBRF and IBRF
quantitatively, we compared the NBRF and IBRF^N^, normalized to the
fixed viewing geometry of NBRF, for the pixel shown by a green circle on
the image of Mongu test site (Figure 4b). The results for the red and
NIR bands are shown in Figure 7 for 2002-2003. The gaps correspond to
the rain seasons with almost continuous cloudiness. The dry season
retrievals show a strong seasonal pattern of vegetation superimposed on
the signals of soil and drying ponds of flood waters. Compared to
IBRF^N^, the response of NBRF to surface change is delayed by as much as
5-8 days. Because cloudiness is very low during the dry season, these
NBRF retrievals were obtained with typically 10-to-16 clear-skies
measurements. Figure 7 also shows a solution obtained with last 4
clear-skies measurements. Obviously, the 4-day solution shows a much
smaller delay (1-2 days) but also a higher noise in a form of several
outliers. With higher cloudiness typical of global performance of ASRVN
and with other uncertainties in the input data including the footprint
variability, the 4-day solution is usually notably noisier than the
current solution which uses all clear-skies measurements in 16-day
Queue. The time delay of the NBRF can be removed by comparing it with
the time series of IBRF^N^, which is very similar between the 4-day
solution and 10-day solution.

As was mention above, variability of footprint with the view angle is
the main source of variability of IBRF. To illustrate this point, Figure
7 shows an IBRF^N^ averaged over an area of 33 pixels (9 km^2^), which
has a substantially lower noise than the 1 km solution. It is
interesting to see that spatial averaging also reduces the temporal bias
between IBRF^N^ and NBRF.

### 4.3 NDVI

The Normalized Differential Vegetation Index (NDVI) is a commonly used
parameter to characterize vegetation canopies. NDVI can be directly
derived from the ASRVN BRF and albedo products. NDVI can be generated
from different reflectance parameters, and its sensitivity to viewing
geometry and atmospheric conditions can change accordingly. Figure 8
shows the NDVI time series for the "green" and "urban" pixels of the
GSFC site as calculated from NBRF, surface albedo, IBRF (representing
different forms of top-of-canopy NDVI) and from TOA reflectance. One can
see that, in agreement with expectations, the variation of NDVI derived
from NBRF and surface albedo is much smaller than that derived from IBRF
and TOA reflectance. The NDVIs derived from the atmospherically
corrected products (NBRF, albedo and IBRF) are also greater than the TOA
NDVI. This is an expected since atmospheric correction tends to reduce
the red band reflectance more than the NIR signal. Comparing the NDVIs
derived from IBRF and NBRF, we find that the NBRF NDVI generally
responds to the seasonal changes with a small delay of several days,
with an exception of cases when larger delay can be traced to a gap in
the AERONET data.

### 4.4 MODIS Terra vs. Aqua

ASRVN creates a sensor-specific time series record of the surface
reflectance over AERONET sites. Figure 9 compares the NBRF time series
for 1.5 years of the "bright pixel" of the Solar Village site. The NBRFs
are very close to each other which suggests a good relative calibration
between the MODIS Terra and Aqua instruments. Generally,
cross-calibration of sensors flying in different orbits is a very
difficult task. With the accumulation of a longer time record and
sufficient global statistics, the ASRVN dataset may become a valuable
source for the cross-calibration analysis of different sensors.

## 5. Concluding Remarks

This paper presented a new operational data collection and processing
system ASRVN, which was initially designed for validation of the surface
reflectance products. ASRVN collects MODIS and MISR L1B data for 5050
km^2^ areas for about 160 AERONET sites, and performs independent
rigorous atmospheric correction using AERONET aerosol and column water
vapor data. The atmospheric correction is achieved by fitting 4-16 -day
(multi-angle) sets of MODIS TOA measurements with theoretical
reflectance accurately parameterized in terms of coefficients of the
LSRT BRF model. The algorithm has a thorough data quality analysis
component, including a cloud mask, an aerosol filter, and the control of
the time series consistency of surface BRF and albedo.

The algorithm is optimized in terms of noise reduction and its ability
to track both seasonal and rapid surface change. During stable surface
conditions, such as periods of maximum greenness during summertime in
the northern hemisphere, a 16-day BRF retrieval gives a solution
characterized by low noise. When the surface changes rapidly (e.g.,
agricultural harvesting), the algorithm gives greater weight to the last
MODIS measurement of the period thus quickly adjusting to the surface
change.

We have also developed a back-up mechanism to cover extended periods of
high cloudiness, when the number of cloud-free MODIS measurements in
16-day Queue is insufficient for the BRF retrieval. In this case, we use
BRF scaling. This approach requires knowledge of the BRF shape which
usually comes from the previous solution stored in the Q-memory. If no
retrievals were made in the past 32 days, the algorithm substitutes the
Q-solution with the historic seasonal BRF, which represents a
quarterly-average BRF over years of retrievals. This method allows us to
rapidly initialize retrievals after long periods of cloudiness and
provide continuous gap-filled imagery of high quality.

The ASRVN suite of products includes three parameters of the LSRT model
(*k^L^*, *k^G^*, *k^V^*), surface spectral albedo, NBRF (a BRF value
computed for a standard viewing geometry, VZA=0^o^, SZA=45^o^), and IBRF
(a BRF value for specific viewing geometry of the last MODIS
measurement). All parameters are produced daily for seven 500 m MODIS
bands at gridded 1 km resolution. We do not store vegetation indices
like NDVI which are easy to produce from the available data.

The ASRVN dataset, including 6 years of MODIS TERRA and 1.5 years of
MODIS AQUA data, is available now as a standard MODIS product (MODASRVN)
which can be accessed through LAADS website (Level 1 and atmosphere
archive and Distribution System,
http://ladsweb.nascom.nasa.gov/data/search.html). The products are
accompanied by a Quality Assurance (QA) flag and color-composite RGB
browse images for the TOA MODIS reflectance, NBRF, *IBRF*, cloud mask
and QA.

The algorithm has high computational efficiency. For example, full ASRVN
re-processing of 6 years of MODIS TERRA data at 100 AERONET sites takes
about 24-30 hours using one processor of a 2.2 GHz workstation.

The results show very stable and reproducible NBRF and NDVI time series
for any given pixel. The main sources of errors in the developed
algorithm are the residual cloudiness and the variation of MODIS pixels
with scan angle, which increase by a factor of eight from the nadir view
to the edge of scan. This effect is partially mitigated by a 1 km
resolution gridding procedure, but it cannot be cancelled entirely. This
source of uncertainty is important in regions with the high surface
heterogeneity at the scale comparable to the grid size of 1 km. In this
regard, expansion of ASRVN with data from the geostationary sensors,
such as the future GOES-R, is expected to produce a higher quality data
set.

The ASRVN applications range from product validation and science
analysis to sensor calibration support and long-term trending and
stability studies. We believe the products from ASRVN fit well into the
CEOS BELMANIP framework and will assist in more reliable and
quantitative intercomparison analysis over the AERONET sites. Recently,
we conducted a validation study of the MISR surface BRF and albedo
products \[11\]. Because ASRVN produces a multi-year record for each
sensor of interest, these data are useful for sensor cross-calibration
analysis \[30\] and detection of long-term calibration trends. The
latter are particularly important for climate applications.

The ASRVN prepares a foundation to support the NPOESS validation program
for the Visible Infrared Imaging Radiometer Suite (VIIRS) which will
replace MODIS and NOAA's Advanced Very High Resolution Radiometer
(AVHRR) as the nation's wide-swath multispectral sensor. The 22-band
VIIRS will provide most of the spectral measurements and capabilities
afforded by MODIS, with bands highly compatible with the existing ASRVN
framework. A constrained pixel growth with the scan angle of the VIIRS
as compared to MODIS is an important factor in the data set quality.
Because NPOESS is an operational program, validation resources will be
significantly more limited than during the EOS era. The NPOESS
Validation Team has defined an initial strategy that strongly emphasizes
the use of operational field network datasets (e.g., AERONET + ASRVN)
over field campaigns. We envision that, similar to MODIS, the
operational processing center will generate VIIRS spatial subsets for
ASRVN in near real time. Several targeted post-launch field campaigns
will be undertaken to verify ASRVN results and document performance.

Our approach can be considered as an indirect validation network for
current MODIS or MISR surface reflectance and associated products (e.g.,
NDVI). If supported by periodic ground measurements over carefully
selected stable homogenous test sites with different surface brightness
level, which would establish an absolute reference for BRF and albedo,
this approach can also be considered a full validation that is easily
expandable to global level given the AERONET global infrastructure.

**Acknowledgements**. The work of Dr. Lyapustin and Dr. Wang was
supported by the NASA EOS Science (Dr. D. Wickland) grant and NASA NPP
(Dr. J. Gleason) grant.

**APPENDIX**

**A-1. Parameterized expression for the TOA radiance**

The algorithm is based on a high accuracy semi-analytical formula
derived with the Green's function method \[25, 31\]. Below, is
atmospheric optical thickness, *S* is spectral extraterrestrial solar
irradiance, and *s*=($\mu = \text{cos}\theta$, ) is a vector of
direction defined by zenith () and azimuthal () angles. The *z*-axis is
pointed downwards, so ~0~\>0 for the solar beam and \<0 for the
reflected beam. The TOA radiance $L(s_{0},s)$ is expressed as a sum of
the atmospheric path radiance (*D*), and surface-reflected radiance
(*L~s~*), directly and diffusely transmitted through the atmosphere:

$L(s_{0},s) = D(s_{0},s)$+$L_{s}(s_{0},s)e^{- \tau|\mu|} + L_{s}^{d}(s_{0},s)$.
(A-1)

The surface-reflected radiance is written as:

$L_{s}(s_{0},s) \simeq$*S*~0~$e^{\frac{- \tau}{\mu_{0}}}${$\rho(s_{0},s)$+${\alpha c}_{0}\rho_{1}(\mu)\rho_{2}(\mu_{0})$} +
$\frac{\alpha}{\pi}\int_{\Omega^{+}}^{}{D_{s}(s_{0},s^{'})\rho(s^{'},s)\mu^{'}ds^{'}}$,
(A-2)

where $D_{s}$ is path radiance incident on the surface, *c*~0~ is
spherical albedo of the atmosphere, and

$\rho_{1}(\mu) = \frac{1}{2\pi}\int_{\Omega^{+}}^{}{\rho(s^{'},s)ds^{'}}$,
$\rho_{2}(\mu_{0}) = \frac{1}{2\pi}\int_{\Omega^{-}}^{}{\rho(s_{0},s)\text{ds}}$.
(A-3)

is a multiple reflection factor, $\alpha = (1 - q(\mu_{0})c_{0})^{- 1}$,
where *q* is surface albedo. The diffusely transmitted surface-reflected
radiance at the TOA is calculated from *L~s~* with the help of 1D
diffuse Green's function of the atmosphere:

$L_{s}^{d}(s_{0},s)$ =
$\int_{\Omega^{-}}^{}{G^{d}(s_{1},s)L_{s}(s_{0},s_{1})\text{ds}_{1}}$.
(A-4)

The function *G^d^* is often called bi-directional upward diffuse
transmittance of the atmosphere. The method of its calculation was
discussed in detail in \[23\]. The surface albedo is defined as a ratio
of reflected and incident radiative fluxes at the surface:

$q(\mu_{0}) = F^{\text{Up}}(\mu_{0}\frac{)}{F^{\text{Down}}}(\mu_{0})$,
(A-5a)

$F^{\text{Down}}(\mu_{0})$=*S*~0~$e^{\frac{- \tau}{\mu_{0}}}$+$\int_{\Omega^{+}}^{}{D_{s}(s_{0},s^{'})\mu^{'}ds^{'}}$=$F_{s}^{\text{Dir}}(\mu_{0}) + F_{s}^{\text{Dif}}(\mu_{0})$,
(A-5b)

$F^{\text{Up}}(\mu_{0}) = {\pi S}_{\lambda}\mu_{0}e^{\frac{- \tau}{\mu_{0}}}q_{2}(\mu_{0}) + \int_{\Omega^{+}}^{}{\mu^{'}q_{2}(\mu^{'})}D_{s}(s_{0},s^{'})ds^{'}$,
$q_{2}(\mu_{0}) = \frac{1}{\pi}\int_{\Omega^{-}}^{}{\rho(s_{0},s)\mu\text{ds}}$.
(A-5c)

These formulas give an explicit expression for the TOA radiance as a
function of surface BRF. The accuracy of the above formulas is high,
usually within a few tenths of a percent \[25\]. Below we will use the
TOA reflectance, which is defined as

$R_{\lambda} = \frac{L_{\lambda}}{(}\mu_{0}S_{\lambda})$. (A-6)

**A-2. Expression for the TOA reflectance using *LSRT* BRF model**

Based on the described semi-analytical solution, we can express TOA
reflectance as an explicit function of parameters of the BRF model. We
are using a semi-empirical Li Sparse -- Ross Thick (LSRT) BRF model
\[32\] as a basis. This is a linear model, represented as a sum of
Lambertian, geometric-optical, and volume scattering components:

$\rho(\mu_{0},\mu,\phi) = k^{L} + k^{G}f_{G}(\mu_{0},\mu,\phi) + k^{V}f_{V}(\mu_{0},\mu,\phi)$.
(A-7)

It uses predefined geometric functions (kernels) $f_{G}$, $f_{V}$ to
describe different angular shapes. The kernels are independent of the
land conditions. The BRF of a pixel is characterized by a combination of
three kernel weights,
$\overrightarrow{K} = \left\{ k^{L},k^{G},k^{V} \right\}^{T}$. The LSRT
model is used in the operational MODIS BRF/albedo algorithm \[9\].

The substitution of equation (A-7) into (A-1 - A-5), and normalization
to the reflectance units gives the following expressions for the
surface-reflected signal (the last two terms of Eq. (A-1)):

$R_{s}(\mu_{0},\mu,\phi)$=
$e^{- \tau\mu_{0}}${$k^{L} + k^{G}f_{G}(\mu_{0},\mu,\phi) + k^{V}f_{V}(\mu_{0},\mu,\phi)$+${\alpha c}_{0}\rho_{1}(\mu)$$\rho_{2}(\mu_{0})$}

\+$\text{αμ}_{0}^{- 1}\left\{ k^{L}E_{0}^{d}(\mu_{0}) + k^{G}D_{G}^{1}(\mu_{0},\mu,\phi) + k^{V}D_{V}^{1}(\mu_{0},\mu,\phi) \right\}$,
(A-8)

$R_{s}^{d}(\mu_{0},\mu,\phi)$=
$e^{- \tau\mu_{0}}${$\lbrack k^{L}G^{\text{av}}(\mu) + k^{G}G_{G}^{1}(\mu_{0},\mu,\phi) + k^{V}G_{V}^{1}(\mu_{0},\mu,\phi)\rbrack$
+

> ${\alpha c}_{0}$$\lbrack k^{L}G^{\text{av}}(\mu) + k^{G}G_{G}^{\text{11}}(\mu) + k^{V}G_{V}^{\text{11}}(\mu)\rbrack$$\rho_{2}(\mu_{0})$}

\+$\text{αμ}_{0}^{- 1}\left\{ k^{L}E_{0}^{d}(\mu_{0})G^{\text{av}}(\mu) + k^{G}H_{G}^{1}(\mu_{0},\mu,\phi) + k^{V}H_{V}^{1}(\mu_{0},\mu,\phi) \right\}$.
(A-9)

The surface albedo is written as:

$q(\mu_{0})$=
$E_{0}^{- 1}(\mu_{0})${~0~$e^{- \tau\mu_{0}}$$q_{2}(\mu_{0}) + k^{L}E_{0}^{d}(\mu_{0}) + k^{G}D_{G}^{3}(\mu_{0}) + k^{V}D_{V}^{3}(\mu_{0})$}.
(A-10)

Different functions of these equations represent different integrals of
the incident path radiance (*D~s~*) and atmospheric Green's function
(*G*) with the BRF kernels. They were described in \[28\] along with the
numerical calculation method. Below, we give only the integral
expressions:

$\rho_{1}(\mu)$=$k^{L} + k^{G}f_{G}^{1}(\mu) + k^{V}f_{V}^{1}(\mu)$,
(A-11)

$\rho_{2}(\mu_{0})$=$k^{L} + k^{G}f_{G}^{2}(\mu_{0}) + k^{V}f_{V}^{2}(\mu_{0})$,
(A-12)

$q_{2}(\mu_{0})$=$k^{L} + k^{G}f_{G}^{3}(\mu_{0}) + k^{V}f_{V}^{3}(\mu_{0})$,
(A-13)

$D_{k}^{1}(\mu_{0},\mu,\phi - \phi_{0})$=$\frac{1}{\pi}\int_{0}^{1}{\mu^{'}d\mu^{'}\int_{0}^{2\pi}{d\phi^{'}}D_{s}(\mu_{0},\mu^{'},\phi^{'} - \phi_{0})f_{k}(\mu^{'},\mu,\phi - \phi^{'})}$,
(A-14)

$D_{k}^{3}(\mu_{0}) = \frac{1}{\pi}\int_{0}^{2\pi}{d\phi^{'}}\int_{0}^{1}{\mu^{'}f_{k}^{3}(\mu^{'})}D_{s}(\mu_{0},\mu^{'};\phi^{'})d\mu^{'}$,
(A-15)

$G^{\text{av}}(\mu)$=$\int_{- 1}^{0}{d\mu}_{1}\int_{0}^{2\pi}{G^{d}(\mu_{1},\mu,\phi - \phi_{1}){d\phi}_{1}}$,
(A-16)

$G_{k}^{\text{11}}(\mu)$=$\int_{- 1}^{0}{f_{k}^{1}(\mu_{1}){d\mu}_{1}}\int_{0}^{2\pi}{G^{d}(\mu_{1},\mu,\phi - \phi_{1}){d\phi}_{1}},$
(A-17)

$G_{k}^{1}(\mu_{0},\mu,\phi - \phi_{0})$=$\int_{- 1}^{0}{d\mu}_{1}\int_{0}^{2\pi}{G^{d}(\mu_{1},\mu,\phi - \phi_{1})f_{k}(\mu_{0},\mu_{1},\phi_{1} - \phi_{0}){d\phi}_{1}}$,
(A-18)

$H_{k}^{1}(\mu_{0},\mu,\phi - \phi_{0})$=$\int_{- 1}^{0}{d\mu}_{1}\int_{0}^{2\pi}{G^{d}(\mu_{1},\mu,\phi - \phi_{1})D_{k}^{1}(\mu_{0},\mu_{1},\phi_{1} - \phi_{0}){d\phi}_{1}}$.
(A-19)

The subscript *k* in the above expressions refers to either
geometric-optical (*G*) or volumetric (*V*) kernels, and the
supplementary functions of the BRF kernels are given by:

$f_{k}^{1}(\mu) = \frac{1}{2\pi}\int_{0}^{1}{d\mu^{'}\int_{0}^{2\pi}{f_{k}(\mu^{'},\mu,\phi^{'} - \phi)}}d\phi^{'}$,
(A-20a)

$f_{k}^{2}(\mu_{0}) = \frac{1}{2\pi}\int_{- 1}^{0}{{d\mu}_{1}\int_{0}^{2\pi}{f_{k}(\mu_{0},\mu_{1},\phi_{1} - \phi_{0})}}{d\phi}_{1}$,
(A-20b)

$f_{k}^{3}(\mu^{'})$=$\frac{1}{\pi}\int_{- 1}^{0}{\mu d\mu\int_{0}^{2\pi}{f_{k}(\mu^{'},\mu,\phi - \phi^{'})}}d\phi$.
(A-20c)

The diffuse and total spectral surface irradiance are calculated from
(A-5b) as:

$E_{0}^{d}(\mu_{0}) = F^{\text{Dif}}(\mu_{0}\frac{)}{(}{\pi S}_{\lambda})$,
$E_{0}(\mu_{0}) = F^{\text{Down}}(\mu_{0}\frac{)}{(}{\pi S}_{\lambda})$.
(A-21)

Let us re-write these equations separating the kernel weights. First,
single-out small terms proportional to the product
$c_{0}\rho_{2}(\mu_{0})$ into the non-linear term:

$R^{\text{nl}}(\mu_{0},\mu)$=${\alpha c}_{0}$$\rho_{2}(\mu_{0})$$e^{- \tau\mu_{0}}${$e^{- \tau|\mu|}$$\rho_{1}(\mu)$ +
$k^{L}G^{\text{av}}(\mu) + k^{G}G_{G}^{\text{11}}(\mu) + k^{V}G_{V}^{\text{11}}(\mu)$}.
(A-22)

Second, collect all remaining multiplicative factors for the kernel
weights:

$F^{L}(\mu_{0},\mu)$ =
($e^{- \tau\mu_{0}}$+$\text{αμ}_{0}^{- 1}E_{0}^{d}(\mu_{0})$)($e^{- \tau|\mu|}$ +
$G^{\text{av}}(\mu)$), (A-23)

$F^{k}(\mu_{0},\mu;\phi)$ =
{$e^{- \tau\mu_{0}}$$f_{k}(\mu_{0},\mu,\phi)$+$\text{αμ}_{0}^{- 1}D_{k}^{1}(\mu_{0},\mu,\phi)$}$e^{- \tau|\mu|}$
+

$e^{- \tau\mu_{0}}$$G_{k}^{1}(\mu_{0},\mu,\phi)$+$\text{αμ}_{0}^{- 1}H_{k}^{1}(\mu_{0},\mu,\phi)$,
*k*=*V*, *G*. (A-24)

With these notations, the TOA reflectance becomes:

$R(\mu_{0},\mu,\phi) = R^{D}(\mu_{0},\mu,\phi)$+$k^{L}F^{L}(\mu_{0},\mu)$+$k^{G}F^{G}(\mu_{0},\mu,\phi)$+$k^{V}F^{V}(\mu_{0},\mu,\phi)$+$R^{\text{nl}}(\mu_{0},\mu)$.
(A-25)

This equation, representing TOA reflectance as an explicit function of
the BRF model parameters, provides the means for an efficient
atmospheric correction.

## References

1.  King, M., & Greenstone, R., *EOS reference handbook: a guide to
    Earth ScienceEnterprise and the Earth Observong System*. Greenbelt,
    MD: EOS Project Science Office, NASA/Goddard Space Flight Center,
    1999, p. 355.

2.  Morisette, J. T., J. L. Privette, C. O. Justice, A framework for the
    validation of MODIS land products. *Remote Sens. Environ*., **83**,
    pp. 77-96, 2002.

3.  Wang, Y., Woodcock, C. E., Buermann, W., Stenberg, P., Voipio, P.,
    Smolander, H., Häme, T., Tian, Y., Hu, J., Knyazikhin, Y. and
    Myneni, R. B., Evaluation of the MODIS LAI algorithm at a coniferous
    forest site in Finland, *Remote Sens. Environ*., 91:114-127, 2004.

4.  Tian, Y., Woodcock, C. E., Wang, Y., Privette, J. L., Shabanov, N.
    V., Zhou, L., Buermann, W., Dong, J., Veikkanen, B., Hame, T.,
    Ozdogan M., Knyazikhin Y., and Myneni, R. B., Multiscale Analysis
    and Validation of the MODIS LAI Product over Maun, Botswana, I.
    uncertainty assessment, *Remote Sens. Environ.*, 83:414-430, 2002.

5.  Tian, Y., Woodcock, C. E., Wang, Y., Privette, J. L., Shabanov, N.
    V., Zhou, L., Buermann, W., Dong, J., Veikkanen, B., Hame, T.,
    Ozdogan M., Knyazikhin Y., and Myneni, R. B., Multiscale Analysis
    and Validation of the MODIS LAI Product over Maun, Botswana, II.
    Sampling strategy, *Remote Sens. Environ*., 83:431-441, 2002.

6.  Garrigues, S, Lacaze R, Baret F, Morisette JT, Weiss M, Nickeson JE,
    Fernandes R, Plummer S., Shabanov NV, Myneni R. B, Knyazikhin Y, and
    Yang W, Validation and Intercomparison of Global Leaf Area Index
    Products Derived From Remote Sensing Data, *J. Geophys. Res*., in
    press.

7.  Baret, F., J. Morisette, R. Fernandes, J.L. Champeaux, R. Myneni, J.
    Chen, S. Plummer, M. Weiss, C. Bacour, G.Derive, , Evaluation of the
    representativeness of networks of sites for the validation and
    inter-comparison of global land biophysical products. Proposition of
    the CEOS-BELMANIP, *IEEE Trans. Geosci. Remote Sens*.,
    44(7)1794-1803, 2006.

8.  Holben, B. N., T. F. Eck, I. Slutsker, D. Tanré, J. P. Buis, A.
    Setzer, E. Vermote, J. A. Reagan, Y. J. Kaufman, T. Nakajima , F.
    Lavenu, I. Jankowiak, A. Smirnov, AERONET-A Federated Instrument
    Network and Data Archive for Aerosol Characterization, *Rem. Sens.
    Environ.*, **66,** 1-16, 1998.

9.  Schaaf, C. B., F. Gao, A. H. Strahler, W. Lucht, X. Li, T.
    Tsang, N. C. Strugnell, X. Zhang, Y. Jin, J.-P. Muller, P. Lewis, M.
    Barnsley, P. Hobson, M. Disney, G. Roberts, M. Dunderdale, C.
    Doll, R. d\'Entremont, B. Hu, S. Liang, and J. L. Privette, First
    Operational BRDF, Albedo and Nadir Reflectance Products from MODIS,
    *Remote Sens. Environ*., 83, 135-148, 2002.

10. Strahler, A. H. and Muller, J. --P., 1999. MODIS BRDF/Albedo
    Product: Algorithm theoretical basis document. NASA EOS-MODIS Doc.,
    Version 5.0, 53 pp.

11. Lyapustin, A., Wang, Y., J. Martonchik, Privette, J., Holben, B.,
    Slutsker, I., Sinyuk, A., Smirnov, A., Local Analysis of MISR
    Surface BRF and Albedo Over GSFC and Mongu AERONET Sites*, IEEE
    Trans. Geosci. Remote Sens*., 44:1707-1718, 2006.

12. Dubovik, O., M. D. King, A flexible inversion algorithm for
    retrieval of aerosol optical properties from Sun and sky radiance
    measurements, *J. Geophys. Res*., **105**, 20,673-20,696, 2000.

13. Dubovik, O., B. Holben, T. F. Eck., A. Smirnov, Y. J. Kaufman, M. D.
    King, D. Tanre, and I. Slutzker, Variability of absorption and
    optical properties of key aerosol types observed in worldwide
    locations, *J. Atmos. Sci*., **59**, 590-608, 2002.

14. Dubovik, O., Sinyuk, A., Lapyonok, T., Holben, B. N., Mishchenko,
    M., Yang, P., Eck, T. F., Volten, H., Munoz, O., Veihelmann, B., van
    der Zande, W. J., Leon, J., Sorokin, M., and Slutsker, I.,
    "Application of spheroid models to account for aerosol particle
    nonsphericity in remote sensing of desert dust", *J* *.Geophys.
    Res*., vol. 111, D11208, doi:10.1029/2005JD006619, 2006.

15. Dubovik, O., A. Smirnov, B. N. Holben, M. D. King, Y. J.
    Kaufman, T. F. Eck, and I. Slutsker: Accuracy assessments of aerosol
    optical properties retrieved from AERONET sun and sky-radiance
    measurements, *J.Geophys. Res*., **105**, 9791-9806, 2000.

16. Lyapustin, A. I., and T. Z. Muldashev, Method of spherical harmonics
    in the radiative transfer problem with non-Lambertian surface, *J.
    Quant. Spectrosc. Radiat. Transfer*, **61**, 545-555, 1999.

17. Lyapustin, A., Interpolation and Profile Correction (IPC) method for
    shortwave radiative transfer in spectral intervals of gaseous
    absorption. *J. Atmos. Sci.*, **60**, 865-871, 2003.

18. Levy, R. C., L. A. Remer, S. Mattoo, E. F. Vermote, Y. J. Kaufman,
    Second-generation operational algorithm: Retrieval of aerosol
    properties over land from inversion of Moderate Resolution Imaging
    Spectroradiometer spectral reflectance. *J. Geophys. Res*., **112**,
    D13211, doi:10.1029/2006JD008711, 2006.

19. Wolfe, R. E., Roy, D. P., and Vermote, E., MODIS Land Data Storage,
    Gridding, and Compositing Methodology: Level 2 Grid. *IEEE Trans.
    Geosci. Remote Sens.*, 36, 1324--1338, 1998.

20. Smirnov A., B.N. Holben, T.F. Eck, O. Dubovik, and I. Slutsker,
    Cloud screening and quality control algorithms for the AERONET data
    base, *Rem.Sens.Env., 73*, 337-349, 2000.

21. Lyapustin, A., Y. Wang, and R. Frey, "An Automatic Cloud Mask
    Algorithm Based on Time Series of MODIS Measurements", *J. Geophys.
    Res*., **113**, D16207, doi:10.1029/2007JD009641, 2007.

22. Li, X. and A. H. Strahler, "Geometric-optical bidirectional
    reflectance modeling of the discrete crown vegetation canopy: Effect
    of crown shape and mutual shadowing", *IEEE Trans. Geosci. Remote
    Sens.*, **30**, 276-292, 1992.

23. Vermote, E.F., A. Vermeulen, 1999, Atmospheric correction algorithm:
    spectral reflectances (MOD09). Algorithm theoretical basis document.
    NASA EOS-MODIS Doc., Version 4.0, 107 pp.

24. Lyapustin, A. I., Atmospheric and geometrical effects on land
    surface albedo, *J. Geophys. Res*., **104**, 4123-4143, 1999.

25. Lyapustin, A., and Yu. Knyazikhin, Green's function method in the
    radiative transfer problem. I: Homogeneous non-Lambertian surface,
    *Appl. Optics*, **40**, 3495-3501, 2001.

26. Lyapustin, A., and Yu. Knyazikhin, Green's function method in the
    radiative transfer problem. II: Spatially heterogeneous anisotropic
    surface, *Appl. Optics*, **41**, 5600-5606, 2002.

27. Tucker, C. J., Red and photographic infrared linear combination for
    monitoring vegetation, *Rem. Sens .Env.*, 8(2): 127-250, 1979.

28. Strugnell, N., and W. Lucht, An algorithm to infer continental-scale
    albedo from AVHRR data, land cover class and field observations of
    typical BRDFs, *J. Climate*, **14**, 1360-1376, 2001.

29. Lucht, W., and P. Lewis, Theoretical noise sensitivity of BRDF and
    albedo retrieval from the EOS-MODIS and MISR sensors with respect to
    angular sampling, *Int. J. Remote Sensing*, **21**, 81-98, 2000.

30. Lyapustin, A., Y. Wang, R. Kahn, J. Xiong, A. Ignatov, R. Wolfe, A.
    Wu, B. Holben, C. Bruegge, Analysis of MODIS-MISR calibration
    differences using surface albedo around AERONET sites and cloud
    reflectance. *Rem. Sens. Environment*, **107**, 12-21, 2007.

31. Lyapustin, A., and Y. Wang, Parameterized Code Sharm-3D for
    Radiative Transfer Over Inhomogeneous Surfaces, *Appl. Optics*,
    **44**, 7602-7610, 2005.

32. Lucht W., Schaaf C. B., Strahler A. H., An algorithm for the
    retrieval of albedo from space using semiempirical BRDF models,
    *IEEE Trans. Geosci. Remote Sens.*, **38**, 977-998, 2000.

Table 1. ASRVN product suite

  ----------- -------------- ---------------- ------------------- ---------------
  **Product   **Data Type**  **Dimensions**   **Descriptions**    **MODIS Product
  Name**                                                          Counterpart**

  CloudMask   DFNT_UINT8     50x50            Cloud mask field,   MOD35 cloud
                                              the definition of   Mask
                                              values is shown in  
                                              Table 2.            

  NBRF        DFNT_FLOAT32   50x50x7          Bi-directional      MOD43B4 Nadir
                                              reflectance factor  BRDF-Adjusted
                                              normalized to       Reflectance
                                              SZA=45^o^ and nadir (NBAR)
                                              view.               

  Albedo      DFNT_FLOAT32   50x50x7          Surface albedo at a A combination
                                              given solar zenith  of MOD43B3
                                              angle in ambient    black-sky and
                                              atmospheric         white-sky
                                              conditions.         albedo weighted
                                                                  with respective
                                                                  relative direct
                                                                  and diffuse
                                                                  incident
                                                                  fluxes.

  IBRF        DFNT_FLOAT32   50x50x7          Instantaneous (or   MOD09 Surface
                                              one-angle) BRF for  reflectance
                                              specific viewing    
                                              geometry of the     
                                              last day of         
                                              observations.       

  Kiso        DFNT_FLOAT32   50x50x7          BRF model           MOD43B1
                                              parameter, the      BRDF/Albedo
                                              isotropic           Model
                                              coefficient.        Parameters

  Kvol        DFNT_FLOAT32   50x50x7          BRF model           MOD43B1
                                              parameter, the      BRDF/Albedo
                                              volumetric          Model
                                              coefficient.        Parameters

  Kgeo        DFNT_FLOAT32   50x50x7          BRF model           MOD43B1
                                              parameter, the      BRDF/Albedo
                                              geometric-optics    Model
                                              coefficient.        Parameters

  QA          DFNT_UINT16    50x50            Quality assurance   N/A
                                              flags. The          
                                              definition of each  
                                              field can be found  
                                              in Table 3.         

  RGB browse  N/A            50x50            The RGB browse      N/A
  image                                       image including TOA 
                                              reflectance, NBRF,  
                                              IBRF, CM and QA.    
  ----------- -------------- ---------------- ------------------- ---------------

Table 2. Cloud mask field definition

  ------------------------------- ---------------------------------------
  **Value**                       **Definition**

  1                               Clear

  2                               Possibly clear

  3                               Possibly cloudy

  4                               Cloudy

  5                               Cloud shadow

  6                               "grey" area around cloud pixel, or
                                  inhomogeneous aerosol.

  10                              Clear water

  11                              Clear snow

  12                              Clear Land (equivalent to 1)

  13                              Clear water (equivalent to 10)
  ------------------------------- ---------------------------------------

Table 3. QA field definition

+----------+------+--------+-----------------------------------------+
| **Bit    | **Bi | **R    | **Bit-code definition**                 |
| field**  | ts** | ange** |                                         |
+----------+------+--------+-----------------------------------------+
| QA       | 0-1  | 00-11  | 00 -- Good quality                      |
| .overall |      |        |                                         |
|          |      |        | 01 -- acceptable quality                |
|          |      |        |                                         |
|          |      |        | 10 -- BRF parameters and NBRF are       |
|          |      |        | filled with previous results, IBRF is   |
|          |      |        | not created                             |
|          |      |        |                                         |
|          |      |        | 11 -- No retrieval is made              |
+----------+------+--------+-----------------------------------------+
| QA.scale | 2    | 0-1    | 0 -- Not scaled                         |
|          |      |        |                                         |
|          |      |        | 1 -- scaled                             |
+----------+------+--------+-----------------------------------------+
| Q        | 3-7  | 0-31   | The number of days since last update of |
| A.nDelay |      |        | Q-memory                                |
+----------+------+--------+-----------------------------------------+
| Q        | 8-9  | 00-11  | 00 -- the BRF consistency parameter     |
| A.status |      |        | *status* =0                             |
|          |      |        |                                         |
|          |      |        | 01 -- *status* =1                       |
|          |      |        |                                         |
|          |      |        | 10 -- *status* =2                       |
|          |      |        |                                         |
|          |      |        | 11 -- *status* \>2 (reliable solution)  |
+----------+------+--------+-----------------------------------------+
| QA.      | 10   | 0-1    | 0 -- no adjacent cloud                  |
| adjCloud |      |        |                                         |
|          |      |        | 1 -- this pixel is adjacent to a cloudy |
|          |      |        | pixel                                   |
+----------+------+--------+-----------------------------------------+
| QA.model | 11   | 0-1    | 0 -- the calculated IBRF is consistent  |
|          |      |        | with the model prediction               |
|          |      |        |                                         |
|          |      |        | 1 -- the calculated IBRF is not         |
|          |      |        | consistent with the model prediction    |
+----------+------+--------+-----------------------------------------+
| QA.cloud | 12   | 0-1    | 0 -- the pixel is clear                 |
|          |      |        |                                         |
|          |      |        | 1 -- the pixel is cloudy                |
+----------+------+--------+-----------------------------------------+

**Figure 1**. Structure of Queue for *ASRVN* processing. The Queue,
designed for the sliding window algorithm, stores up to 16 days of
gridded MODIS observations at 1 km resolution. The CM algorithm uses
MODIS bands 1-7 and bands 31-32, which are stored as Layers
(double-indexed arrays) shown in the upper-left corner. A dedicated
Q-memory is allocated to store the ancillary information for the cloud
mask algorithm in *Refcm* structure, such as a reference clear-skies
image (*refcm*) and results of dynamic Land-Water-Snow classification
(*mask_LWS*). This information is updated with latest measurements (day
L) once given block is found cloud-free, thus adapting to changing
surface conditions. The Q-memory also stores results of previous
reliable BRF retrievals, or atmospheric correction (AC), for MODIS bands
1-7.

![](media/image1.png){width="2.648425196850394in"
height="2.2381889763779528in"}![](media/image2.png){width="2.6570866141732283in"
height="1.3464566929133859in"}![](media/image3.png){width="2.6401574803149606in"
height="1.294488188976378in"}

**Figure 2**. Example of ASRVN processing for GSFC, USA, 2006. The two
left columns show the top-of-atmosphere gridded RGB MODIS TERRA images
for three different sequences of observations. The two left images are
differently normalized to help distinguish clouds. The fourth column
shows the total generated cloud mask (CM legend: clear -- blue, cloud --
red, possibly cloud -- yellow, cloud shadow -- dark red, aerosol filter
-- grey). The last column shows the aerosol optical depth (AOD)
retrieved from the last day of measurement, which is used in aerosol
filter. The AOD is shown using a rainbow palette with pink and red
representing minimum and maximum values, respectively. The third and
fifth columns show the results of ASRVN atmospheric correction (see sec.
3.3).

![](media/image4.wmf){width="0.25078740157480317in"
height="0.16023622047244093in"}

**Figure 3**. Block-diagram of MODIS AC algorithm.

![](media/image5.png){width="1.4826771653543307in"
height="0.5in"}![](media/image6.png){width="1.4826771653543307in"
height="0.5in"}![](media/image7.png){width="1.4826771653543307in"
height="0.5007874015748032in"}![](media/image8.png){width="1.4826771653543307in"
height="0.5in"}![](media/image9.png){width="1.4826771653543307in"
height="0.5in"}![](media/image10.png){width="1.4826771653543307in"
height="0.5in"}![](media/image11.png){width="1.4826771653543307in"
height="0.5in"}![](media/image12.png){width="1.4826771653543307in"
height="0.5007874015748032in"}![](media/image13.png){width="1.4826771653543307in"
height="0.5in"}![](media/image14.png){width="1.4826771653543307in"
height="0.5in"}

\(a\) (b) (c)

**Figure 4**. Temporal dynamics of surface NBRF for a) GSFC in year
2005; b) Mongu site in year 2005 and c) Solar-village site in year 2006.
Shown from left to right are the RGB TOA MODIS gridded images, the RGB
NBRF images and the cloud masks. The true color images are composed from
equally weighted red, green and blue bands. The CM legend is the same as
in Figure 2. The red, yellow and blue dots in panels a) and c) show the
location of selected "green" pixel, "urban" pixel and "bright" pixel,
respectively. The green dot in panel b) shows the location of the test
pixel in Figure 7.

a\)

![](media/image35.wmf){width="5.9944444444444445in"
height="2.6180555555555554in"}

b\)

![](media/image36.wmf){width="5.9944444444444445in"
height="2.6180555555555554in"}

c\)

![](media/image37.wmf){width="5.9944444444444445in"
height="2.6180555555555554in"}

**Figure 5**. NBRF time series for a) green pixel, GSFC site; b) urban
pixel, GSFC site and c) bright pixel, Solar village site. R, G, B color
represents red, green and blue bands, respectively; black is NIR band
and brown is band 7 (2.1m).

![](media/image38.png){width="2.0833333333333335in" height="7.8125in"}

**Figure 6**. Comparison of NBRF and IBRF for GSFC, USA, day 95-122,
2006. Images shown from left to right are: RGB TOA MODIS gridded images,
RGB NBRF images, RBG IBRF images and cloud mask. The CM legend is the
same as in Figure 2.

![](media/image39.wmf){width="7.002777777777778in"
height="3.501388888888889in"}

![](media/image40.wmf){width="7.000694444444444in" height="3.93125in"}

**Figure 7**. Comparison of NBRF (red) and IBRF^N^ normalized to the
standard viewing geometry of NBRF (yellow) for the Mongu site for the
pixel shown by green circle in Figure 4b. Also shown is the NBRF
obtained from 4 last clear-skies days of the Queue (green), and the
IBRF^N^ averaged over an area of 9 km^2^ (blue).

a\)

![](media/image41.wmf){width="5.993055555555555in"
height="2.8402777777777777in"}

b\)

![](media/image42.wmf){width="5.993055555555555in"
height="2.8402777777777777in"}

**Figure 8**. NDVI time series of a) "green pixel" and b) "urban" pixel
for GSFC site.

![](media/image43.wmf){width="5.997222222222222in" height="3.75in"}

**Figure 9**. Comparison of MODIS Terra and Aqua NBRF for a bright pixel
(Solar Village). Here, triangles and solid circles represent Aqua and
Terra NBRF, respectively. Red, green and blue bands are shown with their
respective color, black shows NIR band and brown corresponds to band 7
(2.1m).