## Abstract

This paper advances hyperspectral infrared (IR) radiative transfer techniques for retrieving water (ocean and lake) surface skin temperature from clear-sky radiance observations obtained within the longwave atmospheric window region (800–1000 cm^{−1}). High spectral resolution has optimal potential for multispectral algorithms because of the capability to resolve, and thus avoid, gas absorption lines that otherwise obscure the surface signal in conventional narrowband radiometers. A hyperspectral radiative transfer model (RTM) is developed for varying satellite zenith angles, atmospheric profiles (cloud and aerosol free), surface wind speeds and skin temperatures, with atmospheric column transmittance spectra computed from fast models. Wind speed variations in surface emissivity and quasi-specular reflection are both rigorously accounted for. The RTM is then used for deriving retrieval algorithms based upon statistical and physical methodologies. The statistical method is based upon linear regression analyses of brightness temperatures, whereas the physical method is based upon solution of a linear perturbation form of the IR radiative transfer equation valid for window channels. The physical method is unique in its simplicity: It does not solve for atmospheric profiles, but rather relies upon local linearities about guess transmittances for extrapolating the skin temperature. Both algorithms are tested against independent forward calculations and then used to retrieve water surface skin temperatures from the National Polar-orbiting Operational Environmental Satellite System (NPOESS) Airborne Sounder Testbed-Interferometer (NAST-I) flown on board the NASA ER-2. The results demonstrate the capability of hyperspectral radiative transfer for providing an optimal correction for atmospheric gas absorption (viz., water vapor) from the new suite of environmental satellite IR spectrometers.

## 1. Introduction

Accurate and precise measurements of ocean and lake surface skin temperatures from satellites are required in numerous research and civil applications. Of particular importance are global sea surface temperatures (SST), which are required at a high degree of absolute accuracy (≲0.1 K decade^{−1} system stability) for climate change detection (Allen et al. 1994). Beginning in the early 1970s, retrievals of SST have been performed operationally by the U.S. National Oceanic and Atmospheric Administration's (NOAA's) National Environmental Satellite Data and Information Service (NESDIS) using infrared (IR) detectors flown on board NOAA environmental satellites. Since that time, the accuracy and precision of operational and retrospective products have gradually improved. While SST retrieval is rightly considered a “mature” science, the challenge remains to meet user accuracy requirements through the use of new instruments and retrieval algorithms.

At the present time, there exist numerous instruments and algorithms for estimating surface temperature from space (e.g., Barton 1995). Because the earth's thermal emission peaks in the IR spectrum, and the IR emissivity of water is very close to unity, IR instruments offer maximum signal at the highest spatial resolution.^{1} The usual algorithm approach involves inferential-statistical regression of 2–3 narrowband window-channel brightness temperatures against a known ground truth using data obtained either empirically (e.g., Walton et al. 1998) or from radiative transfer modeling (e.g., Závody et al. 1995; Merchant et al. 1999; Merchant and Harris 1999). These multichannel algorithms provide a correction for molecular water vapor [(H_{2}O)_{υ}] absorption under cloud-free conditions^{2} without requiring atmospheric profile information (McMillin 1975). Among the more notable narrowband instruments are the NOAA Advanced Very High Resolution Radiometer (AVHRR; not designed originally for SST, but provides a data record since 1981; e.g., Walton et al. 1998), the European Remote Sensing Satellites (ERS) Along Track Scanning Radiometer (ATSR/2; designed for SST retrieval with improved detector performance, calibration, and conical-scan for dual-view capability; e.g., Merchant et al. 1999), and the Moderate Resolution Imaging Spectroradiometer (MODIS; offers high spatial resolution with channels optimized for cloud and aerosol characterization). Imagers from geostationary platforms such as the NOAA Geostationary Operational Environmental Satellites (GOES) and the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) Meteosat Second Generation (MSG), have similar spectral channels but offer high temporal resolution.

The above narrowband instruments are each specified to have two spectral channels that divide the 800–1000 cm^{−1} longwave IR (LWIR) window region (commonly referred to as the 11- and 12-*µ*m split window), and an additional window channel in the shortwave IR (SWIR) window region, 2400–2700 cm^{−1}.^{3} The use of the SWIR channel with either of the split-window channels is often referred to as dual-window. It is well documented that retrievals using all three channels are more robust than those using only two (e.g., Hagan 1989; Walton et al. 1998; Merchant and Harris 1999). Unfortunately, however, split-window algorithms are still desirable for daytime retrievals to avoid SWIR solar contamination.^{4} While dual or split-window approaches can yield unbiased results for global data, algorithm-related random and systematic residual errors are observed to occur, particularly in very dry or moist atmospheres, and/or at large satellite zenith angles. Moisture-dependent errors can be significant, especially when considering regions with persistent (H_{2}O)_{υ} loadings. Errors also typically increase with satellite zenith angle as the atmospheric optical pathlength increases and the surface emissivity decreases. This is a crucial problem when considering retrievals from geostationary platforms.

This paper thus reexamines the problem of retrieving accurate and precise water surface skin temperatures, *T*_{s}, with attention given to hyperspectral IR radiative transfer. Physically speaking, there are two essential difficulties to surmount in retrieving *T*_{s}: 1) variable atmospheric attenuation and 2) variable surface radiative properties. Although the first of these includes the effects of opaque and semitransparent clouds, as well as aerosol, this paper will limit its attention to the problem of reducing residual errors associated with gas absorption under globally varying conditions. The problems of cloud detection, or correcting for semitransparent clouds and/or elevated quantities of aerosol (e.g., Merchant et al. 1999; Nalli and Stowe 2002), are acknowledged to be important issues, but they are beyond the scope of this work. High spectral resolution provides the ability to resolve, and thus avoid, gas absorption lines [viz., (H_{2}O)_{υ} and other trace gases] that otherwise obscure the surface signal in conventional narrowband radiometers. These “off-line” hyperspectral channels have become known in the remote sensing community as microwindows. Furthermore, hyperspectral data allow for multiple channels within the traditional split-window region. This capability provides additional degrees of freedom, which may better account for nonlinear variations in water vapor (e.g., Hagan 1989) without solar contamination during daylight. Notable examples of hyperspectral instruments include the Atmospheric IR Sounder (AIRS; on board *Aqua,* launched in 2002), National Polar-orbiting Operational Environmental Satellite System (NPOESS) Preparatory Project Cross-track Infrared Sounder (CrIS) and the MetOp satellite Infrared Atmospheric Sounding Interferometer (IASI; both tentatively planned for 2005), NASA Geostationary Imaging Fourier Transform Spectrometer (GIFTS), and the GOES-R Hyperspectral Environmental Sounder (GOES-R HES).

In section 2, we present a radiative transfer model (RTM) that is used for deriving two different hyperspectral IR retrieval algorithms. The first approach is a conventional inferential-statistical method, whereas the second is a physical retrieval method. The advantages and disadvantages of each approach are considered and discussed. Both retrieval schemes are theoretically validated against forward model calculations in section 4. Finally, in section 5, ocean and lake *T*_{s} are retrieved from radiance spectra obtained by the NPOESS Airborne Sounder Testbed-Interferometer (NAST-I) from the NASA ER-2 during two field campaigns. Retrievals obtained within the vicinity of moored buoys are in reasonable agreement with independent estimates of surface temperature taken from both in situ and polar-orbiting satellite observing platforms.

## 2. IR radiative transfer

The retrieval methods derived in this paper depend upon a realistic, yet practical, radiative transfer model that takes into account the relevant clear-sky physics, including transmittance (Eyre and Woolf 1988), surface emissivity (Wu and Smith 1997), and quasi-specular reflection (Nalli et al. 2001). These model elements are discussed in detail below.

Assuming a plane-parallel, nonscattering, clear-sky atmosphere (i.e., cloud-free with background aerosol levels) with azimuthal symmetry, the IR radiative transfer equation (RTE) for a downlooking detector observing the surface at zenith angle *θ*_{0} and quasi-monochromatic wavenumber *ν* is given by

where *θ* is local zenith incidence angle, *θ*_{0} is the zenith view angle, *R*_{ν}(*θ*_{0}) is the radiance observation (in mW m^{−2} sr^{−1} cm), *ε*_{ν}(*θ*_{0}) is the surface emissivity, *r*_{ν}(*θ, **θ*_{0}) is the bidirectional reflectance (sr^{−1}), *B*_{ν} is the Planck blackbody function, and *T*_{s} is the surface skin temperature. The surface-to-observer path transmittance, 𝒯_{νs}(*θ*_{0}), is defined as

where *τ*_{ν}(*p*_{0}, *p*_{s}) is the vertical optical thickness (or optical depth) of the absorbing species [primarily (H_{2}O)_{υ}, CO_{2}, O_{3}, N_{2}O, CH_{4}, CCl_{4}, and the chlorofluorocarbons CFC11 and CFC12], *p* is pressure altitude, and (*p*_{0}, *p*_{s}) are the observer and surface pressure levels. The upwelling and downwelling atmospheric radiances, *I*^{↑}_{ν}(*θ*_{0}) and *I*^{↓}_{ν}(*θ*), are, respectively, given as

where *T* is the atmospheric layer temperature at pressure altitude *p,* and 𝒯_{ν0} is the observer-to-surface path transmittance (for the *ER-2* flight altitude and higher, 𝒯_{ν0} ≈ 𝒯_{νs} in microwindow channels).

Equations (1)–(4) provide the theoretical basis for IR radiative transfer under cloud-free conditions. Physically, they indicate that there are two fundamental problems to be surmounted in the task of determining *T*_{s} from satellite. The first problem is to model the surface-leaving radiance, prescribed by the term in parentheses on the right-hand side of RTE (1), in a manner that accounts for variations in surface roughness and atmospheric downwelling radiance. Our approach to this problem is discussed in section 2a below. The second problem is to account for the attenuation of this source term through the atmosphere. This is difficult primarily because the profiles of atmospheric water vapor and temperature are highly variable in time and space. This problem is dealt with later in section 3.

### a. Surface-leaving radiance

Treatment of the surface-leaving radiance term is not trivial as is apparent by the complexity of the bidirectional reflectance integral. A method for modeling this term in a manner suitably accurate for SST retrieval is presented here.

#### 1) Surface emissivity

Because the IR emissivity of water is generally ≲0.993 (e.g., Masuda et al. 1988), the radiometric brightness temperature at the surface appears cooler than the actual thermodynamic skin temperature. This is particularly true for sensors operating in spectral window regions, where the deficit is only partially compensated by the reflection of atmospheric radiation (e.g., Smith et al. 1996). A 0.5% departure from a blackbody will yield a significant temperature deficit on the order of 0.3 K or more (e.g., Nalli 1995; Wu and Smith 1996).

In recent years, remote sensing scientists have sought to model and account accurately for the nonunity emissivity of the ocean surface using numerical wave-facet models (e.g., Masuda et al. 1988; Závody et al. 1995; Watts et al. 1996; Wu and Smith 1997; Nalli and Smith 1998; Merchant et al. 1999). In this work, lookup tables (LUT) of Wu and Smith (1997) emissivity spectra are used: *ε*_{ν} = *f*(*ν, **θ*_{0}, *V*), where *V* is the mean surface wind speed and *θ*_{0} is the detector view angle. For lake surfaces, the refractive indices of pure water are used, whereas for ocean surfaces a salinity correction is applied (see Nalli et al. 2001).

#### 2) Surface-reflected radiance

Modeling the surface-reflected radiance presents a more challenging problem due to the dependence on the hemispherical distribution of downwelling atmospheric radiation (e.g., Watts et al. 1996; Nalli et al. 2001). In this work we employ a quasi-specular model developed by Nalli et al. (2001) for calculating the reflected radiance in window channels. For convenience in retrievals, the full hemispherical model calculation is greatly reduced by defining a reflection-diffusivity angle, *θ*_{ν}, such that

where *R*_{νs}(*θ*_{0}) is radiance reflected from the surface, and

The diffusivity angle is simply the mean angle of downwelling radiance reflected from the surface into the detector (dependent on both surface roughness and atmospheric opacity). To allow a convenient determination of reflected radiance for a given atmosphere and sea state, LUT are computed for *θ*_{ν} by finding the zeros of Eq. (6) for a range of wavenumbers, viewing angles, wind speeds, and uplooking slant-path transmittances, 𝒯_{ν0}(*θ*_{0}), hence *θ*_{ν} = *f*[*ν, **θ*_{0}, *V*, 𝒯_{ν0}(*θ*_{0})].

### b. RTM for retrievals

The combined use of LUT for *R*_{νs} and *ε*_{ν} provides an efficient and suitably accurate specification of the surface radiometric properties as functions of surface wind speed and atmospheric transmittance (Nalli et al. 2001). Substituting Eq. (6) and the model emissivity *ε*_{ν}(*θ*_{0}, *V*) into Eq. (1), and including spectral random instrument noise, leaves a simplified retrieval model RTE

where *N*_{ν} is a Gaussian random number and *σ*_{ν} is the instrument noise equivalent radiance (NER), calculated as the standard deviation of a sample of ambient-temperature, calibration-blackbody spectral deviations from the Planck function. Although the NAST-I NER typically exhibits spectral variations, a nominal constant value of 0.2 mW m^{−2} sr^{−1} cm is assumed in the simulations for simplicity. From Eq. (8), we can proceed to derive multispectral retrieval algorithms for *T*_{s}.

## 3. Multispectral retrieval methods

In this section, statistical and physical multispectral retrieval methodologies are developed based upon RTM Eq. (8). The hyperspectral window channels used for both the physical and statistical retrievals are first described below.

### a. Hyperspectral window channels

#### 1) Microwindows

Microwindow channels are selected according to empirically determined threshold values of absorbing-gas component surface-to-observer (50 hPa) University of Maryland, Baltimore County (UMBC) fast model transmittances for H_{2}O continuum (Clough 1995), H_{2}O lines, O_{3}, and “fixed gases” (viz., CO_{2} and CFCs). Channels are kept only if all their component transmittances are greater than the threshold values. Figure 1 shows the component transmittances calculated for the sample atmosphere and the consequent selection of microwindow channels. Depending on the threshold selection, the number of microwindows selected by this process can easily total *n*_{ν} > 100 in the 800–1000 cm^{−1} region. Note that the hyperspectral microwindows selected from this process captures the smoothly varying spectral signature of the H_{2}O continuum absorption across the window region.

#### 2) Miniwindows

Although there are *n*_{ν} > 100 channels available to us, much of the information content is redundant and, indeed, not necessary for retrieving *T*_{s}. Redundancy, however, is still desirable because of the level of spectral random noise typically present in IR spectrometer data. Noise obscures the information content in microwindow channels and can lead to instabilities in retrieval algorithms. One practical way of minimizing the effects of noise is to average together adjacent clusters of microwindow channels (A. Harris 2003, personal communication) to yield what will be referred to as “miniwindows.” While doing this, however, it is still important to restrict the miniwindows to narrow spectral bandpasses so as not to invalidate quasi-monochromatic radiative transfer. As of this writing, we find that *n*_{ν} = 9 miniwindows works well, although there is nothing special about this number. Note that miniwindows defined in this manner are much “cleaner” than typical narrowband filter channels and do not require any specification of a spectral response function.

### b. Statistical retrieval

Statistical algorithms are based upon an observed relationship between channel brightness temperatures (predictors) and surface temperature. Although it is possible to explain physically how such observed relationships can occur using a truncated radiative transfer equation (e.g., McMillin 1975), the discovery (e.g., Anding and Kauth 1970) and rigorous justification of any such relationship is rooted in statistical analysis. In practice, there are two approaches for deriving statistical equations: 1) empirically against coincident in situ buoy matchups (e.g., McClain et al. 1985; Walton et al. 1998), or 2) synthetically from RTM calculations based on training profiles (e.g., Merchant et al. 1999). Apart from simplicity, the primary advantage of conventional statistical-empirical approaches is the ability to account implicitly for, in a mean sense, all potential sources of uncertainty, including surface and atmospheric effects, as well as instrument noise and calibration. The disadvantages, however, include the difficulty in acquiring high-quality, globally representative buoy matchup samples (e.g., Nalli and Stowe 2002), along with the inherent limitations of linear regression equations for explaining all observable variance. Radiative transfer based statistical retrievals effectively circumvent the first of these disadvantages, but must still rely upon a parametric linear equation derived from an accurate RTM as is done below.

A multispectral statistical retrieval for miniwindow channels is given by

where *T*_{i} is the brightness temperature for miniwindow channel *i.* The coefficients *b*_{i}(*θ*_{0}) are determined from statistical regression of *T*_{S} against *T*_{i}, derived from model RTE (8) and a global training data sample (described below), for a discreet set of zenith angles spanning the scanning range of the instrument. Retrievals are then obtained at intermediate *θ*_{0} by cubic spline interpolation of *b*_{i}. Note that Eq. (9) is identical to that used for the ATSR instrument (e.g., Merchant et al. 1999), except it is applied to multiple hyperspectral observations without multilook capability.

As stated earlier, statistical approaches work well under limited ranges of atmospheric conditions, but random and systematic errors are known to occur under extreme water vapor conditions due to global nonlinearities (e.g., Hagan and Nalli 2001). Multispectral retrievals could therefore benefit from a priori information about the atmospheric condition (e.g., Steyn-Ross et al. 1993; Barton 1995). Indeed, this is the motivation behind approaches such as the nonlinear SST (e.g., Walton et al. 1998), which employs an empirical nonlinear predictor, or alternatively by using “regional” or moisture-dependent coefficients in Eq. (9) (e.g., Závody et al. 1995). In this work, however, we propose a physical approach that accounts explicitly for atmospheric profile variations by seeking a direct solution for *T*_{s} based upon a simplified form of RTE (8).

### c. Physical retrieval

To solve RTE (8) for *T*_{s}, we first turn our attention to the surface emission term (first term in parentheses on right side), which consists of three potential unknowns factored together: *ε*_{ν}(*θ*_{0}), *B*_{ν}(*T*_{s}), 𝒯_{νs}(*θ*_{0}). Separating these requires a priori information about at least one of the variables. Over land surfaces, emissivity can vary substantially and is poorly known (spatially, temporally, and spectrally), while a priori atmospheric data are readily available. In these circumstances it has been desirable to either separate emissivity and temperature given accurate a priori atmospheric information (e.g., Kahle and Alley 1992; Xie 1993; Gu and Gillespie 2000), or to devise a simultaneous solution for both (e.g., Ma et al. 2002). However, for the ocean the opposite is true: Emissivity is more uniform and predictable, but atmospheric data are typically sparse. Given the LUTs described in section 2a(1), we will assume the emissivity and reflected radiance spectra to be known quantities. Guess estimates of surface wind speed and atmospheric transmittance necessary for the LUTs can be obtained from various sources and will be elaborated on more later.

#### 1) Atmospheric correction

Theoretically speaking, were atmospheric profiles of temperature and water vapor well known a priori, transmittance and radiance profiles can be calculated, and RTE (8) could be solved for *T*_{s} (the only remaining unknown) from a single spectral observation. Such an approach would constitute a single-channel physical retrieval. However, this method cannot provide a reliable solution because of strong dependence on the accuracy of atmospheric profiles available in practice (Nalli and Smith 1998; Hagan and Nalli 2001). Therefore, drawing from the original work of Smith (1970) and McMillin (1975), a simple multispectral solution for *T*_{s} is derived in the following paragraphs.

First, applying Eqs. (4) and (3) to (8) (excluding noise) and expressing the radiances at a mean wavenumber *v* (e.g., Smith et al. 1972), yields

where *T*_{B} is the observed channel brightness temperature given by the inverse Planck function, *T*_{B}(*ν, **θ*_{0}) = *B*^{−1}_{ν}[*R*_{ν}(*θ*_{0})]. Nalli (2000) demonstrated quantitatively that the errors introduced by “normalizing” radiances to *v* in RTE (10) are negligible when applied to a small spectral interval (e.g., 800–1000 cm^{−1}). This was attributed to the effective cancellation of the cumulative spectral radiance departures from *v*.

An iterative forward model (after Smith 1970) for the upwelling radiance is expressed as

where the radiance quantities for the *n*th iteration are denoted with the superscript (*n*). Assuming we can obtain guess estimates of the atmospheric profiles and *T*_{s} [discussed below in section 3a(3)], the channel transmittances are assumed to be known approximately and the initial state (*n* = 0) can be computed. This is subtracted from a single iteration (*n* = 1) to yield a perturbation form of the RTE:

where

It is necessary to interject here that, strictly speaking, the accuracy of the guess transmittances directly depend upon the accuracy of the guess water vapor profile, as well as the transmittance model. Thus, for any given channel, the validity of Eq. (12) is expressly contingent upon suitable guess water vapor profile information. However, as will be demonstrated below, a one-iteration calculation involving multiple channels does not require perfectly known profiles to achieve an accurate and precise solution for surface temperature.

Proceeding then, we apply the mean value theorem (after McMillin 1975) to Eq. (12) over the vertical gradient of transmittance for the upwelling and downwelling radiance streams

which leads to

where *T*_{a↓}, *T*_{a↑} are mean lower-tropospheric temperatures (over the channel weighting functions) for up and downlooking views, respectively. Substituting Eqs. (18) and (19) into (10) results in the following equation:

Disregarding the small channel-dependent variations of *δB*_{v}(*T*_{a↑}) and *δB*_{v}(*T*_{a↓}), the resulting system of equations can be solved for the three unknowns, namely *δB*_{v}(*T*_{s}), *δB*_{v}(*T*_{a↑}), and *δB*_{v}(*T*_{a↓}) given *n* ≥ 3 spectral window channels. However, because of the small inherent magnitude of the reflection term [the second term in Eq. (20)], it is more practical to assume that the downwelling atmospheric radiance perturbation is approximately equal to that for the upwelling (after Nalli and Smith 1998), that is

This approximation is most accurate for transparent microwindow channels and when temperature profile deviations vary slowly with altitude {i.e., *δB*_{v}[*T*(*p*)] ≈ constant}.

where the “known” parameters *x*_{ν} are given by

We note that Eqs. (22)–(23) are similar to the Nalli and Smith (1998) equations (17), (19), and (27), with the exception that the forward radiances *B*_{v}[*T*^{(0)}_{B}(*ν, **θ*_{0})], and “known” parameters *x*_{ν}, contain physics more appropriate for larger zenith view angles (i.e., downwelling transmittances calculated at the reflection-diffusivity angle *θ*_{ν} instead of the specular angle *θ*_{0}).

#### 2) Retrieval algorithm

Equation (22) may be expressed in matrix form as

where *j* = {1, 2}, *δ*y_{ν} denotes an *n* × 1 column vector of observed minus modeled radiances, 𝗫_{νj} denotes an *n* × 2 matrix of “known” parameters [these being *ε*_{ν}(*θ*_{0}, *V*), 𝒯_{νs}(*θ*_{0}), and 𝒯_{νs}(*θ*_{ν})], and *δ*b_{j} denotes an unknown 2 × 1 column vector of deviations from the guess state for the surface and atmospheric emissions. The matrix 𝗫_{νj} is derived for each retrieval location given a guess surface wind speed [to compute *ε*_{ν}(*θ*_{0}, *V*)], and guess profiles [to compute 𝒯_{νs}(*θ*_{0}) and 𝒯_{ν0}(*θ*_{ν})]. Given *n* > 2 window observations, the overdetermined linear system (24) is readily solved for *δ*b_{j} by the method of least squares.^{5} From Eq. (14), the retrieved water surface skin temperature is taken to be

where *δB*_{v}(*T*_{s}) ≡ *δb*_{1} is determined from the solution to (24).

Retrieval equations (22)–(24) explicitly model the radiative processes indicated by RTE (8) for spectral window regions where the surface signal is most evident. The solution of (24) is not to be confused with standard methods of deriving regression coefficients for an assumed parametric statistical model (as was done in section 3b). The methodology is unique in this sense and is properly classified as a physical retrieval algorithm (see Kidder and Vonder Haar 1995). Furthermore, by restricting ourselves to window channels, it is possible to solve (24) for surface temperature with far less computational effort than would a full simultaneous iterative solution of the RTE. Although high spectral resolution is advantageous, only two spectral measurements are mathematically required for solving (24). Hence, theoretically speaking, the algorithm could be used with conventional narrowband channel radiometers, including the AVHRR, ATSR, and MODIS (see Nalli and Smith 1998; Hagan and Nalli 2001, for studies that demonstrate this capability).

Physically speaking, the retrieval may be conceptualized as an extrapolation of radiance perturbations to zero optical depth (i.e., unit transmittance), as the solution *δB*_{v}(*T*_{s}) is simply the *y* intercept of Eq. (22) (see Nalli and Smith 1998, for an illustration of this concept). Provided that the spectral-relative characterization of the channel transmittances is well defined (i.e., the transmittance model is accurate), absolute errors in transmittance spectra (due to erroneous guess water vapor profiles) do not severely impact the extrapolation derived from the multispectral linear Eq. (24). This assertion is demonstrated in section 4 with forward modeling calculations.

The guess profiles and *T*^{(0)}_{s} used in the physical retrieval can be obtained from a variety of sources. The following section suggests one possible approach whereby window channel radiance observations are themselves used to retrieve guess profiles.

#### 3) Guess profile retrievals

Guess temperature and water vapor profiles can be retrieved from window channel observations using a simple statistical method (e.g., Smith and Woolf 1976). The equation for a given atmospheric pressure level *p*_{l} and satellite zenith angle *θ*_{0} is expressed as (B. Huang 1999, personal communication)

where *X* denotes either temperature or water vapor, *n* is the number of spectral channels, an overbar denotes the mean, and *a*(*ν*_{i}, *p*_{l}, *θ*_{0}) are coefficients to be determined from model calculations of *T*_{B}(*ν*_{i}, *θ*_{0}) using a representative sample of training profiles (described in section 4). Equation (26) is expressed in matrix form as

where the solution for coefficients 𝗔 given by

are determined for a set of discreet zenith angles spanning the scanning range of the instrument. Profile retrievals are then obtained at intermediate angles by cubic spline interpolation of 𝗔.

It might at first seem doubtful that profiles retrieved in this manner using only window channels could be useful given weighting functions that are limited solely to the lower troposphere. However, it is precisely these weighting functions that are germane to the radiative transfer within the window channels. In this sense, the retrieved lowest layer air temperature can also be used as the guess skin temperature, *T*^{(0)}_{s}, necessary in retrieval equation (22). The potential for obtaining more accurate guess profiles in practice, either from satellite radiance measurements across the CO_{2} and H_{2}O absorption bands, or from a numerical weather prediction model, are beyond the scope of this paper.

## 4. RTM calculations

To determine the theoretical accuracy of the algorithms, high-resolution microwindow radiance spectra are calculated using fast transmittance models appropriate for NAST-I. The standard downlooking fast model used for simulation was developed at the University of Wisconsin (UW)—Madison (after Eyre and Woolf 1988). For downwelling calculations, an uplooking model was specifically developed courtesy of H. Woolf, UW—Madison.

### a. Data samples

A quality-controlled, marine sample (*N* = 9024) of atmospheric profiles is extracted from a global 1998 profile dataset (H. Woolf, UW—Madison, 1999, personal communication) to represent conditions of the global ocean. These profiles are used in Eq. (8) to simulate NAST-I radiance observations, complete with spectral random noise, for various satellite zenith angles, surface wind speeds, and surface temperatures.

A subsample of *n* = 873 profiles is used for training statistical equations (27) and (9), with the training *T*_{s} for (27) taken as three deviations, Δ*T* = {−2, 0, +2} K, from the lowest layer air temperature. Radiance spectra are calculated for four surface wind speeds *V* = {0, 5, 10, 15} m s^{−1} at five satellite zenith angles, *θ*_{0} = {0°, 33.56°, 44.42°, 51.32°, 56.25°}.^{6} Thus, four wind speeds and 3 *T*_{s} deviations yields 4 × 3 × 873 = 10 476 forward radiance calculations for deriving statistical retrieval coefficients at each of the five *θ*_{0}.

An independent profile data sample (*n* = 301), along with uniquely different values of *T*_{s}, *V*, and *θ*_{0}, is used for validation of the retrieval methods. For these independent forward calculations, the radiance “observations” are simulated based upon RTE (1) with the “true” bidirectional surface reflection explicitly calculated [i.e., see Nalli et al. 2001, Eqs. (24) and (A5)] without using approximation (6). These calculations include spectral random noise as in Eq. (8). The “true” *T*_{s} are modeled as a random number deviation (*µ* = 0; *σ* = 1.5 K) from the lowest layer air temperature. The “true” surface wind speeds and zenith angles are taken to be *V* = {2.5, 7.5, 12.5} m s^{−1} and *θ*_{0} = {24.62°, 45.43°, 55.15°},^{7} respectively. This yields 3 × 301 = 903 forward calculations for each *θ*_{0}.

### b. Retrieval error analyses

#### 1) Guess profile error

It is first desirable to estimate the physical retrieval's sensitivity to guess profile error. This is achieved by simulating profile errors as perturbations (constant with altitude) from the true profiles taken from the independent data sample. Erroneous guess temperature profiles, *T*′(*p*), are modeled as *T*′(*p*) ≡ *T*(*p*) + *δT*, where *T*(*p*) is the true profile and *δT* = {−8, −7, … , +7, +8} K. Erroneous guess (H_{2}O)_{υ} mixing ratio profiles, *r*^{′}_{w}(*p*), are modeled with fractional perturbations *δr*_{w}/*r*_{w}(*p*) = {−0.5, −0.4, … , +0.4, + 0.5} g kg^{−1}, that is, *r*^{′}_{w}(*p*) ≡ *r*_{w}(*p*) + *δr*_{w}. The corresponding systematic and random errors of the *T*_{s} retrievals for the 301 independent profiles are then computed. A robust estimate of the systematic error (bias) is given by the median, whereas the random error (scatter) is given by the robust standard deviation (RSD), here defined as the median absolute deviation divided by the factor 0.6745. The reason for using these descriptive statistical measures (as opposed to the mean and standard deviation) is that they provide robust estimates of central tendency and scatter that minimize sensitivity to outliers (e.g., Merchant and Harris 1999; Nalli and Stowe 2002). The results for wind speed *V* = 2.5 m s^{−1} are given in Fig. 2 (row-wise for each zenith angle), with bias and random errors given in the left- and right-hand columns, respectively. Overall, the physical retrieval shows relatively little sensitivity to guess profile errors in Fig. 2. Negatively biased *r*_{w}(*p*) tend to yield negative *T*^{(1)}_{s} bias with some scatter, whereas positively biased *r*_{w}(*p*) yield lesser positive bias with greater scatter. *T*(*p*) errors have little to no impact on *T*^{(1)}_{s} bias, but can exacerbate random error when coupled with erroneous *r*_{w}(*p*).

In the error analysis that follows [section 4b(2)], the guess profiles necessary for the physical retrieval are retrieved from microwindow channels^{8} and Eq. (27). Figure 3 shows the corresponding guess profile error statistics for wind speed *V* = 2.5 m s^{−1}. In general, the guess profile retrievals do not exhibit any strong dependencies with respect to *V* (figures not shown here) or *θ*_{0}. Profile errors tend to increase with altitude as would be expected from microwindow weighting functions concentrated in lower troposphere.

Recall that the physical retrieval also requires an a priori estimate of surface wind speed, *V*. In the retrieval simulations that follow, errors in guess *V* are modeled as a Gaussian random error with a standard deviation of 2 m s^{−1} (the accuracy threshold of wind retrievals from passive MW instruments; N.-Yu Wang 2002, personal communication).

#### 2) Skin temperature error

Simulated *T*_{s} retrievals for the independent data sample are obtained from physical retrieval equations (24) and (25), and independently from statistical retrieval equation (9). Table 1 summarizes the results according to both conventional and robust statistical measures. Statistical significance for the robust measures are inferred from the Efron bootstrap technique (Efron 1979), as shown in Fig. 4. The bootstrap technique generates large numbers of artificial samples through random resampling (with replacement) of an observed sample. This provides a significance measure for sample statistics without making any assumptions about the population distribution. The independent data retrievals are randomly resampled 1000 times to generate histograms of retrieval bias (navy–blue–cyan) along with scatter about the mean (yellow–orange–maroon). Three colors are used for each histogram to distinguish the results obtained under the three wind speed regimes. From Fig. 4, the following conclusions may be drawn about the global performance of the physical and statistical algorithms.

The retrievals from both algorithms show little to no global bias resulting from the error sources considered in this work.

^{9}The statistical retrieval exhibits a small dependence on surface wind speed at the larger two zenith angles. Wind speed dependence is largely mitigated with the physical algorithm.

The physical algorithm yields less random error than the statistical algorithm, most notably for larger zenith angles. While these differences may appear small, they amount to an improvement in precision by approximately 30%. Furthermore, these are statistically significant results as demonstrated by comparing the yellow-orange-maroon histograms obtained from each approach and noting the small overlap.

To examine retrieval dependencies on atmospheric (H_{2}O)_{υ} content, Fig. 5 shows scatterplots of retrieval error, *δT*_{s}, versus total column precipitable water Σ_{H2O} for all three model wind speeds. In general, scatter is seen to increase from lower to higher water vapor and from lower to higher satellite zenith angles. A close inspection reveals that the physical algorithm is generally more robust (i.e., exhibits less localized bias and scatter) with respect to global variations in water vapor. We also expect that residual errors could be further reduced given more accurate guess profiles (cf., Figs. 2 and 3) The statistical retrievals, on the other hand, exhibit nonlinearity and become significantly biased negative for moist atmospheres (Σ_{H2O} ≳ 40 mm). Modeling studies conducted by Hagan and Nalli (2001) for narrowband split-window instruments demonstrate even more pronounced scatter and nonlinearity. These model results confirm that LWIR-based linear regression equations do not adequately model the variance observed over the global range of (H_{2}O)_{υ}.

## 5. Experimental case studies

Field observations of high-resolution IR radiance spectra are derived from NAST-I (Cousins and Smith 1997), an experimental, scanning Fourier transform spectrometer flown on board the NASA ER-2 research aircraft. The NAST-I obtains spectral measurements in three bands ranging over 645–2700 cm^{−1} with an unapodized resolution of ≃0.25 cm^{−1}. The surface “footprint” at nadir is ≃2.6 km from the *ER-2* flight altitude of 20 km. The cross-track mirror scan of ±45° yields a surface observational swath width of ≃46 km. The onboard calibration is accurate to within 1%, with the system random noise equivalent temperature (NEdT) ≃0.3 K in the LWIR band. Note that this system noise is significantly higher than what was ideally assumed in the retrieval simulations from the previous section. This additional random error should be kept in mind when considering the NAST-I retrievals in the following paragraphs.

NAST-I data used for analysis were obtained over lake and ocean during two field campaigns in 1999. The first of these was the Winter Experiment (WINTEX), which was conducted out of Madison, Wisconsin, from 15 March to 2 April 1999, under winter–spring midlatitude climate conditions. The second experiment was conducted out of the NASA Flight Research Center, Wallops Island, Virginia, under midlatitude summer conditions during August 1999. Both experiments included clear-sky overpasses of National Data Buoy Center (NDBC) moored buoys, which provided an independent, hourly in situ measurement of bulk water surface temperature at ≃0.5–1.0-m depth. NDBC moored buoys are specified to be accurate to within ±1°C (Hamilton 1986; Meindl and Hamilton 1992). This measurement uncertainty is in addition to any water temperature gradients (vertical and horizontal) as well as temporal undersampling. Therefore, we supplement the buoy measurement with remotely sensed AVHRR Pathfinder retrievals extracted from the AVHRR Pathfinder matchup database (PFMDB; Kilpatrick et al. 2001). These retrievals are closely matched in space and time with the NDBC buoys. The AVHRR Pathfinder retrieval algorithm is a refined statistical method derived retrospectively using empirical regression based upon the PFMDB. Consequently, it is important to keep in mind that the AVHRR retrievals are statistically “tuned” to the buoys.

During WINTEX the aircraft data were obtained over Lake Michigan, whereas data from Wallops were collected over the Atlantic Ocean off the coast of Long Island, New York. Guess profiles used in the *T*_{s} physical retrievals were obtained using Eq. (26) with coefficients derived from the global 1998 training dataset described in section 4a. Unlike simulated data, retrieval accuracy using experimental data will necessarily depend upon the accuracy of the transmittance model. Thus, for these case studies we opted to use the downlooking NAST-I fast model developed by UMBC because it includes CFC absorption bands whereas the UW—Madison model does not. In spite of this, it was still found that 2 miniwindow channels within the 935–965 cm^{−1} region degraded all retrievals (both *T*_{s} and guess profiles) relative to the coincident in situ data. We speculate that this may have resulted from the inherent difficulty of modeling the complex absorption lines (CO_{2} and H_{2}O) in this region (see Fig. 1, “H_{2}OL” and “Dry” components). All retrievals were improved after omitting these two channels.

To demonstrate the derived atmospheric correction from the physical algorithm under moist atmospheric conditions, Fig. 6 shows the spatial distributions of NAST-I observations taken during the Wallops experiment before and after the retrieval. Mapped in Fig. 6a are the spectral mean brightness temperatures for microwindow channels. Limb darkening is evident in this plot, resulting from the increased atmospheric pathlength, as well as the decreased surface emissivity, that occurs at larger zenith view angles. Also, the mean microwindow temperatures are cooler than the reported buoy temperature by approximately 2–3 K. By contrast, the retrieved skin SSTs in Fig. 6b do not exhibit any noticeable limb darkening, and the absolute temperatures are in much closer agreement with the buoy temperature. Significant mesoscale features are also apparent in the retrieved SSTs, with large horizontal gradients of approximately 2–3 K over approximately 20 km present. Such gradients are known to occur and probably result from subsurface temperature structures associated with the dynamics of the region. Furthermore, the Gulf Stream can cause the formation of mesoscale eddies, which also contribute to spatial variability. Note that there has been no spatial smoothing or analysis applied to the NAST-I temperature maps.

In a similar manner, lake *T*_{s} from WINTEX are shown in Fig. 7 that include three separate overpasses of an operational NDBC buoy in Lake Michigan. Evident again are mesoscale features, comparable in size to those seen in the Atlantic Ocean during Wallops. A persistent localized warm spot is also apparent in the vicinity of the buoy. This feature is distinct throughout the three overpasses over the time period 0.75–2.3 UTC, and is reasonably consistent at different scan angles.

Table 2 summarizes NAST-I results from both field campaigns, with reference given to concurrent buoy observations as well as AVHRR Pathfinder retrievals extracted from the PFMDB. Although the AVHRR Pathfinder retrievals extracted from the PFMDB. Although the AVHRR Pathfinder buoy matchups have tight space–time constraints (i.e., within ±30 min and ±0.1° lat–lon of buoy observations; Kilpatrick et al. 2001), it is not possible to obtain perfect overlap between the buoy, NAST-I, and AVHRR observations. Thus, we have sought the best possible matchup in space and time between the Pathfinder AVHRR and NAST-I. Because NAST-I has a higher spatial resolution than AVHRR global area coverage (GAC) data, a spatial mean is calculated for NAST-I field of view (FOV) falling within a (4 km)^{2} extraction box (roughly corresponding to a GAC pixel) centered on the AVHRR Pathfinder retrieval. For the Wallops overpass, only one NAST-I FOV fell within the extraction box, whereas for WINTEX, the number of FOV totalled 3, 2, and 2, for each successive flight leg, respectively. The first column of Table 2 indicates the detector-buoy matchup under consideration (i.e., NAST-I or AVHRR). The second and third columns show the date and time periods for each of the matchups. The fourth and fifth columns show the remotely sensed *T*_{s} from physical and statistical algorithms, *T*^{(1)}_{s} and *T*_{S}, respectively. The sixth column shows the observed channel-mean brightness temperatures, *T*_{B}. The final three columns show the buoy observations for bulk water surface temperature, *T*_{w}, air surface temperature, *T*_{a}, and mean surface wind speed, *V*, respectively.

From the table, it would appear that the water surface temperature retrievals from both NAST-I and AVHRR are biased positive relative to the buoys. For WINTEX, the NAST-I physical retrieval exhibits nearly the same degree of bias as the AVHRR Pathfinder retrieval (≃+0.8 K), whereas the NAST-I statistical retrieval exhibits a more pronounced bias (≃+0.9 to +1.1 K). The larger bias in the statistical retrieval may be due to regression coefficients overcorrecting for the dry continental atmosphere. For Wallops, the NAST-I physical and statistical retrievals are in very close agreement with the buoy (+0.06 and +0.08 K, respectively), whereas the AVHRR Pathfinder is biased warmer (+0.21 K).

Because both NAST-I and AVHRR independently exhibit warm biases, we speculate that there may exist skin–buoy temperature differences during the overpass times. A positive skin–buoy difference is possible when considering the time of day and the meteorological data reported by the buoys. Minnett et al. (1984) have pointed out the potential difficulty for satellite validations against buoys during the late afternoon arising from the development of a diurnal thermocline. During conditions of high insolation and low wind speed, a stratified surface layer can develop yielding a large vertical temperature gradient that is highly variable in space and time. During Wallops, given the time of day [≃1900 UTC = 1400 local standard time (LST)], moderate wind speed (*V* ≈ 4 m s^{−1}), and negligible air–sea temperature difference (*T*_{a} − *T*_{w} = 0 K), it is quite possible that a warm surface layer caused the skin temperature to be slightly warmer than the temperature at 0.6-m depth. For WINTEX, although the time of day was after sunset (≃0100–0200 UTC = 1900–2000 LST), the presence of significant positive air–sea temperature differences (1.3–1.7 K), along with low surface wind speeds (3.0–3.5 m s^{−1}), may also have led to an apparent warm skin temperature relative to the buoy. These skin-buoy uncertainties underscore the need for accurate and precise radiometric *T*_{s} as ground truth for validation [e.g., the Marine Atmospheric Emitted Radiance Interferometer (M-AERI), see Smith et al. 1996; Minnett et al. 2001). Other potential sources of error include space–time matchup variabilities, erroneous guess profile estimates (for the physical retrieval), buoy measurement errors, and instrument noise, calibration and RTM errors. More extensive field validation is clearly needed using the new suite of satellite high-resolution spectrometers (e.g., AIRS, CrIS, IASI, GIFTS, and GOES-R HES).

## 6. Conclusions and future work

It is now three decades since the original pioneering work on multispectral satellite remote sensing of sea surface temperature (e.g., Anding and Kauth 1970; Smith et al. 1970; McMillin 1975). Since that time, notable advances have occurred with innovations in the engineering of remote sensors and in the science of radiative transfer modeling. Advances in instrumentation often lead to corresponding developments in radiative transfer models and retrievals, which then in turn provide steering for future instrument specifications. Although the basic multispectral approach for correcting clear-sky (cloud and aerosol free), IR window channel brightness temperatures for atmospheric water vapor is well known, we reassert that the global problem of ocean and lake surface temperature retrieval is not yet solved to the specification of user accuracy thresholds.

Toward this end, this work has illustrated a distinct advantage of hyperspectral data, namely the unique capability to select *n* > 100 spectral channels (i.e., microwindows) within the 800–1000 cm^{−1} LWIR window region that avoid the spectral absorption bands of IR optically active gases [viz., (H_{2}O)_{υ}, CO_{2}, O_{3}, CFCs]. A radiative transfer model was developed for microwindow channels and used to derive statistical and physical retrieval algorithms. The RTM accounted rigorously for surface emissivity and quasi-specular reflection of downwelling atmospheric radiance. Atmospheric column transmittance spectra were efficiently computed from fast transmittance models. The physical algorithm, based upon a linear perturbation form of the RTM, required only rough estimates of atmospheric temperature, water vapor and surface wind speed. Retrieval simulations showed that although both hyperspectral algorithms were theoretically robust to variations in column water vapor (using the same microwindow channels), the physical algorithm minimized moisture-dependent errors.

Ocean and lake *T*_{s} were then retrieved from the NAST-I instrument flown onboard the NASA *ER-2* during two field experiments during 1999. Retrievals within the vicinity of NDBC moored buoys in the Atlantic Ocean and Lake Superior compared reasonably both with in situ buoy measurements as well as retrievals extracted from the AVHRR Pathfinder matchup database. More extensive validation of the algorithms will be necessary using satellite hyperspectral instruments (e.g., AIRS, CrIS, IASI, GIFTS, and/or GOES-R HES).

The inherent simplicity of the proposed algorithms would facilitate operational implementation within NOAA/NESDIS orbital processing systems. The modeling results demonstrate the capability of high spectral resolution for water vapor correction under globally diverse profile conditions. We acknowledge, however, that in practice the retrieval accuracies will also be contingent upon the accuracy of the transmittance model, instrument characterization, and cloud and aerosol detection algorithms. Aerosols may be accounted for, either physically or statistically, by incorporating ancillary information about the aerosol optical depth into the retrieval equations (e.g., Nalli and Stowe 2002), by selecting channels that can infer a correction, or by synergistically employing passive MW data. Development of aerosol correction methodologies for hyperspectral data will be the subject of ongoing research.

## Acknowledgments

This research was supported under NASA Grants NAS5-31375 (H. E. Revercomb, UW—Madison), NAS5-31377, the Wisconsin Space Grant Consortium, and the NOAA/NESDIS Ocean Remote Sensing Program. We are indebted to H. Woolf (UW—Madison) for his development of the uplooking NAST-I fast-model and compilation of the global 1998 profile data. We are also grateful to A. Harris (University of Maryland, College Park) and B. Huang, R. Knuteson, B. Howell, P. van Delst, C. Sisko, D. Tobin, and D. DeSlover (all UW—Madison) for their expert technical advice and/or assistance. J. Vazquez (JPL/Caltech) helped us in resolving questions pertaining to PFMDB data available on the Internet. Finally, we extend our appreciation to the three anonymous reviewers who provided many excellent comments and insights that have enhanced the paper's scope and quality.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

Current affiliation: NOAA/NESDIS/ORA, Washington, D.C

*Corresponding author address:* Dr. Nicholas R. Nalli, Cooperative Institute for Research in the Atmosphere, NOAA/NESDIS E/RA3, 5200 Auth Rd., Room 603-9, Camp Springs, MD 20746-4304. Email: Nick.Nalli@noaa.gov

^{1}

Although passive microwave (MW) based retrievals are receiving more attention recently due to the insensitivity of microwaves to clouds, aerosols, and water vapor, the focus and discussion in this paper is restricted to IR-based retrievals. Because of the distinct advantages of both IR and MW, a synergistic instrument approach may eventually be desirable for an optimal, blended SST product.

^{2}

Because clouds are typically opaque (i.e., zero transmittance) in the IR (with the exception of certain high-level ice clouds), all IR retrieval algorithms must necessarily assume some degree of cloud-free radiative transfer, even in partly cloudy situations. Therefore, the problem of cloud detection, a difficult undertaking in its own right, is usually treated as a problem distinct from retrievals as is done in this work.

^{3}

The GOES M–Q satellites do not have split-window capability, and MODIS has 3 SWIR channels for SST.

^{4}

A daytime, dual-window SST algorithm for *GOES-12* has been devised by screening regions of sunglint and then accounting for diffuse solar contributions (C. Merchant 2002, personal communication).

^{5}

Note that an exact solution can be obtained given *n* = 2 spectral observations.

^{6}

The zenith angles are chosen as five equidistant points in sec(*θ*_{0}) bounded by 1 and 1.8.

^{7}

Similar to the training data, the zenith angles are chosen as three equidistant points in sec(*θ*_{0}) bounded by 1.1 and 1.75.

^{8}

Guess profiles can also be retrieved from miniwindows, but significant water vapor profile biases were found to result for profiles with large deviations from the mean profile. Because the instrument noise is perfectly characterized in simulation, profile retrievals are more accurate when using microwindows. As mentioned elsewhere in this paper, it may be necessary in practice to obtain more reliable guess profiles using a more rigorous method than the simplistic approach employed here.

^{9}

We recognize, however, that for radiative transfer based retrievals, errors in the instrument calibration and/or transmittance model, as well as residual clouds and aerosol, will probably result in global biases if not properly accounted for in practice (topics beyond the scope of this paper).