Atmospheric Measurement Techniques A sea surface reflectance model for ( A ) ATSR , and application to aerosol retrievals

A model of the sea surface bidirectional reflectance distribution function (BRDF) is presented for the visible and near-IR channels (over the spectral range 550 nm to 1.6 μm) of the dual-viewing Along-Track Scanning Radiometers (ATSRs). The intended application is as part of the Oxford-RAL Aerosols and Clouds (ORAC) retrieval scheme. The model accounts for contributions to the observed reflectance from whitecaps, sun-glint and underlight. Uncertainties in the parametrisations used in the BRDF model are propagated through into the forward model and retrieved state. The new BRDF model offers improved coverage over previous methods, as retrievals are possible into the sunglint region, through the ATSR dual-viewing system. The new model has been applied in the ORAC aerosol retrieval algorithm to process Advanced ATSR (AATSR) data from September 2004 over the south-eastern Pacific. The assumed error budget is shown to be generally appropriate, meaning the retrieved states are consistent with the measurements and a priori assumptions. The resulting field of aerosol optical depth (AOD) is compared with colocated MODIS-Terra observations, AERONET observations at Tahiti, and cruises over the oceanic region. MODIS and AATSR show similar spatial distributions of AOD, although MODIS reports values which are larger and more variable. It is suggested that assumptions in the MODIS aerosol retrieval algorithm may lead to a positive bias in MODIS AOD of order 0.01 at 550 nm over ocean regions where the wind speed is high. Correspondence to: A. M. Sayer (sayer@atm.ox.ac.uk)


Introduction
The Intergovernmental Panel for Climate Change (IPCC) has identified aerosols as among the most uncertain contributions to radiative forcing (Penner et al., 2001, Forster et al., 2007).As approximately 70% of the Earth's surface is covered by water, the accurate determination of aerosol loadings over ocean is critical to assess direct and indirect aerosol radiative effects.In the visible and near-infrared (nIR) spectral domains the ocean surface is dark, particularly compared to typical land surfaces, meaning the proportional atmospheric contribution to the signal measured by imaging radiometers at the top-of-atmosphere (TOA) is higher for the same aerosol loading.However, typical oceanic aerosol loadings are low (see, for example, Smirnov et al., 2009), meaning the surface contribution is non-negligible.An exception to the rule of the ocean being dark is found in sun-glint, whereby solar and satellite geometries lead to regions where the surface is very bright, typically in the tropics for nearnadir-viewing instruments.Parametrisations of sun-glint are largely based on the approach of Cox and Munk (1954a), and most aerosol retrieval algorithms use a glint formulation to identify and mask out glint-affected regions before processing.This has the effect of reducing the spatial coverage of the derived aerosol dataset, particularly in the tropics.
A notable exception to this is given by O' Brien and Mitchell (1988), who relied on the predictable spatial variation of surface reflectance within large cloud-free portions of the sun-glint region to peform aerosol and wind speed retrievals from Advanced Very High Resolution Radiometer (AVHRR) data.This methodology has, however, seen little application since.
Multiangle imaging instruments such as the Along-Track Scanning Radiometers (ATSRs), Multiangle Imaging Spec-troRadiometer (MISR) and POLarization and Directionality Published by Copernicus Publications on behalf of the European Geosciences Union.
of the Earth's Reflectances (POLDER) allow for an improved representation of surface anisotropy in aerosol retrieval algorithms, although over oceans measurements from a single viewing geometry have been considered sufficient to derive useful aerosol information.The treatment of surface reflectance in some of these (single-view or multiview) algorithms is described below.Typically, the primary quantity retrieved is the aerosol optical depth (AOD) at a midvisible wavelength; many algorithms use a fixed surface reflectance, and some a fixed aerosol type.Regional differences exist in ocean aerosol climatologies from different sensors (for example, Thomas et al., 2010).Validation of aerosol optical depth over the open ocean is difficult; many land and coastal regions are well-represented by ground-based measurements taken by the AErosol RObotic NETwork (AERONET, Holben et al., 1998), while the Maritime Aerosol Network (MAN, Smirnov et al., 2009) of shipborne open ocean measurements is spatially and temporally sparser.It is likely that some of the differences between these satellite climatologies arise from the assumptions made about the ocean surface reflectance.Partially as a result of the comparative darkness of the ocean as compared to land surface reflectances, the various algorithms as summarised below tend to show little change since their early versions.
The over-ocean aerosol retrieval algorithm for the MODerate resolution Imaging Spectroradiometer (MODIS) is described first by Tanré et al. (1997), and the same basic algorithm is applied in the current Collection 5 of the dataset (Remer et al., 2005).The methodology is adapted from Koepke (1984), who defined glint, whitecap (foam) and underlight (scattering from dissolved pigments) contributions to the surface reflectance.The glint formulation of Cox and Munk (1954a) is used with a fixed wind speed, whitecap and underlight contribution, and a glint threshold is defined in which no retrievals are performed (unless heavy dust loading is detected, in which case the retrieval is attempted).The overall surface reflectance is fixed in the algorithm.Sediment masks are used to remove pixels of high sediment loading, which are not accounted for by the reflectance algorithm.The MISR ocean aerosol retrieval algorithm (Martonchik et al., 1998) also uses the method of Koepke (1984).
Veefkind and de Leeuw (1998) use a similar algorithm and fixed surface reflectance to retrieve AOD over ocean as a mixture of two aerosol types (anthropogenic and maritime) from ATSR-2 data.The nadir and forward views are used independently, and a comparison between the two views can be used as a consistency check.The study noted that errors in the TOA reflectance arising from an incorrect wind speed could lead to errors of 0.04-0.16 in nadir-view-derived AOD.As typical open ocean optical depths may be of this order (Smirnov et al., 2009), this is a significant possible error source.More recently, this algorithm has been applied by Bennouna et al. (2009) to Spinning Enhanced Visible and Infra-Red Imager (SEVIRI) data.A previous version of the algorithm described here and applied to ATSR-2 data, involv-ing a nadir-view aerosol retrieval algorithm, is detailed by Thomas et al. (2009b) and Thomas et al. (2010).This took a similar approach for sea surface reflectance although allowed the absolute magnitude to vary (while fixing the spectral shape of the surface).Two-channel aerosol retrieval algorithms are presented for AVHRR by Higurashi and Nakajima (1999) and Mishchenko et al. (1999): these use fixed glint-based surface reflectances calculated as described in Nakajima and Tanaka (1983) and Mishchenko and Travis (1997), respectively.They also consider the impacts of the simple reflectance model on the retrieved AOD and Ångstrom exponent, noting that it can be significant for cases of high wind or pigment concentrations, or low aerosol loadings.
Dedicated ocean colour sensors such as the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) focus on the retrieval of parameters such as chlorophyll-a concentration and treat aerosol as part of an atmospheric correction term.Approximations made in surface reflectance models, such as the "black pixel approximation" (that the water-leaving radiance in the nIR is negligible), have been shown to negatively impact upon the quality of retrieved ocean colour parameters (Siegel et al., 2000).Sano (2004) describes an algorithm for retrieval of AOD, Ångstrom exponent and aerosol refractive index from POLDER reflectance and polarisation measurements at 670 nm and 865 nm, making use of the black pixel approximation along with the glint formulation of Cox and Munk (1954a).
This work describes a new algorithm for the calculation of sea surface reflectance, drawing upon the methodology of Koepke (1984).The intended application is as part of the Oxford-RAL Aerosol and Clouds (ORAC) scheme.This is discussed here in the context of aerosol retrievals, although the model may also be applied in the case of opticallythin cloud (where the surface contribution at TOA is nonnegligible).As each sensor is different, previous assumptions must be reevaluated, and more recent work taken into account, to create a model suitable for the ATSRs.By modelling accurately the contributions from different sources, retrievals are possible within the sun-glint region, which increases the possible coverage of aerosol retrievals as compared to existing algorithms.Additionally, potential biases in the retrieved aerosol fields arising from neglect of accounting for foam and underlight are avoided, and all four of the visible/nIR channels on the ATSRs may be used.Furthermore, in ORAC, unlike many previously-described algorithms, the surface reflectance is not fixed, adding some flexibility in those cases in which the assumed surface reflectance is incorrect.
2 Overview of the ORAC-(A)ATSR aerosol retrieval

The (A)ATSR instruments
The ATSR series consists of three instruments: ATSR-1 (aboard ERS-1), launched in 1991, ATSR-2 (aboard ERS-2), launched in 1995, and the Advanced ATSR (AATSR, aboard Envisat), launched in 2002.The ATSRs were primarily designed for measurement of sea surface temperature (Závody et al., 1995).While ATSR-1 measured radiance at one wavelength in the near-infrared and three in the thermal infrared part of the spectrum, ATSR-2 and AATSR have an additional three channels in the visible region.It is these visible channels which are key to the instruments' ability to provide data suitable for aerosol retrievals, and so ATSR-2 and AATSR, referred to from here as (A)ATSR, are considered here.
ERS-2 and Envisat are in Sun-synchronous polar orbit with a mean local solar equatorial crossing time of 10:30 a.m.(ERS-2) or 10:00 a.m.(Envisat) for the descending node.The ATSR instruments are unique in that they use two views (near-simultaneous in time) with differing path lengths to discriminate between radiance from the surface and radiance from the atmosphere.(A)ATSR measures at seven channels in the visible and infrared; at present the first four (centred near 550 nm, 660 nm, 870 nm and 1.6 µm, known as channels 1-4, respectively) are used in the aerosol retrieval scheme.The additional bands are centred near 3.7 µm, 11µm and 12 µm.
The shortwave quantity reported by (A)ATSR for a given channel is an approximation of the spectral bidirectional reflection factor, the Sun-normalised radiance, R TOA (θ s ,φ s ;θ v ,φ v ), which is defined where θ s ,φ s denote the illumination (solar) zenith and azimuth angles and θ v ,φ v the corresponding angles of view (the sensor), respectively.A channel is defined between wavelengths λ 1 ,λ 2 to have response (λ).Finally L r λ is the radiance measured by the instrument and E i λ is the TOA downward solar irradiance.
The area sampled by (A)ATSR consists of two curved swathes: a nadir view, looking down at zenith angles from 0 • -22 • , and a forward view inclined between 53 • -55 • to the normal to the surface.There are 555 pixels across the nadir swath (with a size of about 1 km 2 at the centre) and 371 across the forward swath (with a size of about 1.5 km 2 at the centre).During each scan cycle the satellite moves approximately 1 km onward with respect to the Earth's surface; after around 150 s the satellite has moved such that nadir view samples the same region, giving two views of the scene with differing path lengths.Global coverage is achieved every 3-6 days depending on location.ATSR-2 operates in a narrow-swath mode over much of the ocean, reducing coverage by approximately half, due to data-downlinking restrictions from the ERS-2 platform.
(A)ATSR has an on-board visible calibration system consisting of an opal diffuser which views the Sun once per orbit.This, together with vicarious calibration against stable bright ground targets, means that the visible channel reflectances are known to an accuracy of 2-3% (Smith et al., 2002(Smith et al., , 2008)).

The ORAC retrieval
ORAC is an optimal estimation (OE) retrieval (Rodgers, 2000) making use of Levenburg-Marquardt iteration (Levenberg, 1944;Marquardt, 1963) to find the most probable state of the surface and atmosphere given measurements and a priori information.The measurement vector consists of the TOA reflectances for the nadir and forward views of the first four channels.The retrieved state parameters (the "state vector") are the aerosol optical depth at 550 nm (τ 550 ), the aerosol effective radius (r e ) and the surface bihemispherical reflectance at each of the four channels used (R dd,1 , R dd,2 , R dd,3 and R dd,4 ).The AOD is reported at 550 nm as this is the commonly-used standard; the derived AOD may, however, be referenced to any wavelength, and is obtained from all measurements simultaneously.
For computational speed, cloud-free forward and nadirview data are typically averaged to a 10 km sinusoidal grid before the ORAC retrieval is performed.This averaging to a coarser resolution is known as "superpixelling".From here, the term "ground scene" is taken to refer to the data, superpixelled or not, used for an individual retrieval.However, ORAC can in principle be performed at any resolution.
The robust statistical basis of OE provides the following advantages: 1. Estimates of the quality of the retrieval solution (the retrieval "cost") for each ground scene.This is essentially an error-weighted χ 2 test of the fit to the measurements at the retrieval solution, which provides a level of confidence as to the results of any one retrieval.
2. Estimates of the random error on each retrieved parameter for each ground scene.These arise through knowledge of the uncertainty on the measurements and any a priori data, propagated through the forward model.
3. The ability (but not requirement) to use any a priori data available on the state parameters.The model described in this work provides an a priori for the surface bihemispherical reflectance.
The retrieval forward model, presented in Thomas et al. (2009a), calculates the TOA reflectance for a given viewing geometry and state vector.It makes use of precalculated lookup tables (LUTs) of atmospheric transmission and reflectance using the DISORT radiative transfer code (Stamnes et al., 1988).A selection of aerosol models are used in the www.atmos-meas-tech.net/3/813/2010/Atmos.Meas.Tech., 3, 813-838, 2010 retrieval, corresponding to typical continental, desert, maritime or urban aerosol, using aerosol components drawn from the OPAC database of Hess et al. (1998), and additionally a model for biomass burning aerosol drawn from Dubovik et al. (2002).These models consist of mixtures of aerosol components, and different effective radii are obtained by altering their mixing ratios during the retrieval.Generally, the retrieval is attempted for each aerosol type.The most likely aerosol type may be chosen either by considering the model which resulted in the best fit to the measurements (the lowest retrieval cost), or using other available information external to the retrieval.-Strub et al. (2006) noted that in remote sensing terms relating to reflectance were often misunderstood or applied ambiguously or incorrectly.They defined nomenclature for nine types of reflectance, using the framework of Nicodemus et al. (1977), corresponding to incoming and outgoing radiation that is either directional, conical or hemispherical.

Schaepman
The relevant geometric notation used throughout this work is given in Table 1.For clarity and conciseness of notation, spectral variability of the reflectances is implicit in the definitions and so omitted in the notation.
The most fundamental quantity is the bidirectional reflectance distribution function (BRDF), denoted in the ORAC retrieval by R bb : This defines the BRDF in terms of the proportion L r of the incident irradiance E i reflected from direction (θ s ,φ s ) into direction (θ v ,φ v ).In this case, the point of incidence is the Sun and point of reflection is the satellite sensor.It has units of sr −1 , and as a ratio of infinitesimal quantities it (and other directional reflectances) may not be directly observed.In general use the term is defined as a surface property, although a TOA BRDF could also be defined as a conceptual analogue to R TOA .The BRDF is integrable over angles to obtain the other reflectance quantities given by Schaepman-Strub et al. (2006).Unitless reflectance factors are defined as the ratio of observed radiant flux to the radiant flux reflected under the same geometric conditions by an ideal Lambertian surface, such that the bidirectional reflectance factor (BRF) is the BRDF multiplied by π (Schaepman-Strub et al., 2006).
The closest observable equivalent to the BRF is the biconical (or conical-conical) reflectance factor (BCRF or CCRF), obtained by integrating the BRF over solid angles ω to generate cones of incident and reflected light.Conical quantities become a good approximation for the related directional qualities when the solid angles of the cones are small.In this case the solid angle subtended by the Sun is small, as is the instrument's instantaneous field of view (IFOV) of 1/777 rad.For ambient lighting conditions there will be atmospheric contributions from diffusely-scattered light and absorption.These effects lead to the need for an "atmospheric correction" for ground-based sensing applications, or conversely they provide the "signal" for atmospheric sounding; i.e. the biconical reflectance observed at the TOA may have significantly different spectral and angular characteristics to the biconical reflectance just above the surface.Through optimal estimation, ORAC extracts the information about both from the TOA measurements.
To account for the mixture of direct and diffuse illumination the ORAC retrieval forward model (Thomas et al. , 2009a) treats the direct and diffuse contributions to TOA reflectance with separate terms, subjecting them to different reflectances at ground, and different atmospheric transmittances.The surface reflectances required for direct and diffuse radiance may be derived from the BRDF.Hence it becomes necessary to define three types of surface reflectance in the forward model: 1.The surface BRDF, R bb .This describes the reflection of the direct solar beam into the viewing angle, and is a function of both solar and viewing angles.The BRDF is different for each of AATSR's viewing geometries.This is assumed equivalent to the CCRDF and so no integration over solid angle is performed.
2. The directional-hemispherical reflectance (DHR), R bd .This describes the diffuse reflection of the direct beam over the whole hemisphere (or alternatively direct reflection of incoming diffuse radiance), and is a function of the solar angle.The short time delay between the forward and nadir views means that the solar angle and hence DHR are effectively identical for both views.This is sometimes referred to as black-sky albedo, as incoming illumination comes from a sole direction.
3. The bihemispherical reflectance (BHR), R dd .This describes the reflection of diffuse downwelling radiation, assumed isotropic.Hence it is independent of the geometry, and is the quantity retrieved by the retrieval algorithm.This is sometimes referred to as white-sky albedo, as illumination arises from the whole of the sky.
In this notation the subscript b indicates a direct beam reflectance and d a diffuse reflectance; the DHR R bd , for example, denotes an incoming direct beam being diffusely reflected.Given an analytical description of R bb , the DHR for a given solar zenith angle may be obtained by integration over all satellite viewing zenith and relative (solar-satellite) azimuth angles: This may then be integrated over all solar zenith angles to obtain the BHR: For the ocean surface, Gaussian quadrature integration with 4 points in each angular dimension is sufficient to obtain the DHR and BHR to 3 significant figures from the BRDF (the glint contribution is precalculated with a higher number of points, as discussed later).The BHR at each wavelength used are retrieved by the ORAC scheme, but there is insufficient information to also retrieve the full BRDF from the measurements.Therefore BRDF models are used to generate R bb , R bd and the a priori R dd .The ratios R bb : R dd and R bb : R bd are fixed in the aerosol retrieval, such that when R dd is scaled in an iterative step in the retrieval then these ratios are used to scale R bb and R bd by the corresponding factor.This work describes the sea surface BRDF as calculated in ORAC.

The three components of R bb
The model described in this work draws on the heritage of Koepke (1984).An implementation of the Koepke (1984) description of surface reflectance, focusing on the 400 nm-700 nm spectral range, is in the 6S radiative transfer code described by Vermote et al. (1997).Koepke (1984) describes R bb as being composed of three terms representing different sources of upwelling irradiance.Firstly, light can be reflected off whitecaps in the rough ocean surface; secondly, it can be  reflected off the foam-free portion of the surface.The contributions from these two factors will depend on the roughness of the sea surface, which is determined by the wind speed.Thirdly, light penetrating the surface can be scattered back up into the atmosphere by molecules within the body of water.The combination of these terms leads to the relationship where f wc ρ wc is the contribution to reflectance from whitecaps; ρ gl represents the sun glint; and ρ ul denotes the "underlight" term from radiance reflected just below the surface of the water.This is represented schematically in Fig. 1.Although these components represent reflectances, they are denoted using ρ instead of R for clarity.These three components are dealt with individually due to their differing directional and spectral variability, as summarised in Table 2.The terms ρ gl and ρ ul are weighted against by a factor of (1-f wc ), where f wc is the fractional cover of whitecaps, as specular reflectance and underlight are taken to arise from only the foam-free portion of the surface.The formulation for ρ ul includes a correction to account for light lost due to glint reflection at the surface (see Sect. 6).
Whitecaps are where the ocean appears bright due to the action of wind creating a foam.The simplest of the three components of R bb , their only dependence is on wind speed and wavelength.The contribution of whitecaps to reflectance is the product of the proportion of the surface covered by whitecaps (f wc ) and their average reflectance (ρ wc ).Koepke (1984) treated whitecap reflectance in the visible region as constant with wavelength, although noted that in the nearinfrared it might be expected to decrease due to absorption by water molecules.More recent coastal (Frouin et al., 1996) and open ocean (Nicolas et al., 2001) work suggests a reflectance of about 0.4 at shorter wavelengths, decreasing by about 40% at 850 nm and 85% at 1.65 µm.These ratios have been adopted here for use at the nearby (A)ATSR channels, with reflectance at 550 nm and 660 nm assumed equal to 0.4.Kokhanovsky (2004) develops a physical model for whitecap reflectance, which is then parametrised in terms of the (spectral) water absorption and a spectrally neutral coefficient.This coefficient is determined on a case-by-case basis from several measurement sets, including Frouin et al. (1996).This model suggests that the whitecap reflectance may vary globally.The adoption of global values based on Frouin et al. (1996) and Nicolas et al. (2001) introduces in most cases negligible error into the calculation of R bb as the whitecap fraction is generally low, and the variability among the experimental cases studied in Kokhanovsky (2004) is small compared to uncertainties in the whitecap reflectance (Frouin et al., 1996;Nicolas et al., 2001) and whitecap fraction.
The whitecap fraction is here parameterised in terms of wind speed, w, by a simple power law according to the method of Monahan and Muircheartaigh (1980).The fractional cover of whitecaps is given by with the caveat that f wc cannot be greater than 1.It should be noted that determination of f wc is complicated and various formulations based on wind speed and other environmental factors have been developed.An overview of some of these methods is given by Anguelova and Webster (2006).
The method of Monahan and Muircheartaigh (1980) is used as it has been widely adopted (such as Koepke, 1984) and requires only easily-available wind speed data.In the ORAC retrieval scheme, 6-hourly 10 m winds at 1 degree resolution from the European Centre for Medium-range Weather Forecasts (ECMWF), linearly interpolated in space and time, are used throughout.
The contribution of whitecap reflectance to R bb as a function of wind speed is shown in Fig. 2. As it lacks geometric dependence, the contribution of whitecaps to R bb , R bd and R dd are the same.The global mean wind speed for 2004, sampled at AATSR overpass times, is shown in Fig. 3.For wind speeds of approximately 10 ms −1 and higher f wc nm, the green line the contribution at 870 nm, and the blue line the contribution at 1.6 µm.

39
Fig. 2. Whitecap coverage and contribution to the BRDF as a function of wind speed.The dashed black line indicates the whitecap coverage f wc .The red line shows the contribution f wc ρ wc to R bb at 550 nm and 660 nm, the green line the contribution at 870 nm, and the blue line the contribution at 1.6 µm.and f wc ρ wc are considerable (10 −3 −10 −2 , except at 1.6 µm).Such high wind speeds are found polewards of 45 • , with typical ocean wind speeds elsewhere in the range 5-8 ms −1 , corresponding to ρ wc around 10 −4 −10 −3 .There are several sources of uncertainty with this section of the algorithm: -There is a large uncertainty of up to 50% in the spectral reflectance of whitecaps (Frouin et al., 1996;Nicolas et al., 2001;Kokhanovsky, 2004).
- Anguelova and Webster (2006) reveal that different parameterisations of f wc can lead to estimates differing by up to an order of magnitude, as a simple dependence on wind speed is inadequate to explain observed variability.This may be a significant source of error in high-wind or low-chlorophyll environments away from the sun-glint region.

Glint reflectance
The contribution ρ gl results from rays of light striking the wind-ruffled sea surface and being specularly reflected in the observer's direction.It is calculated using the model of Cox and Munk (1954a) and Cox and Munk (1954b), from which key equations are reproduced here for completeness.More complicated than whitecaps, glint depends strongly on geometry, wind speed and wind direction, and weakly on wavelength.

Calculation
The algorithm defines a coordinate system (P , X, Y, Z) such that P is the observed point on the surface and Z the altitude with P Y in the direction of the Sun and P X in the direction perpendicular to the Sun's plane.The surface slope is defined by the two components where θ s < π/2 and θ v > 0. In reality, the slope distribution will be anisotropic and dependent on wind direction χ w .The axes are rotated clockwise from the north by χ w to define a new coordinate system (P ,X ,Y ,Z) where P Y is parallel to the wind direction, and where the slope components may be re-expressed: Following Cox and Munk (1954a), the probability distribution of surface facets p(Z x ,Z y ) is required to calculate the glint reflectance.The original work provided coefficients for 3 parametrisations for p: -Wind-isotropic (dependent only on absolute wind speed).
-Wind-anisotropic, with an additional Gram-Charlier series correction term.
Recently, Zhang and Wang (2009) evaluated these parametrisations, along with other work drawing on the heritage of Cox andMunk (1954a) (specifically Wu, 1972;Mermelstein et al., 1994;Shaw and Churnside, 1997;Shifrin, 2001;Ebuchi andKizu, 2002, andBréon andHenriot, 2006) using MODIS measurements.It was found that the anisotropic model (without the Gram-Charlier series) of Cox and Munk (1954a), and the model of Bréon and Henriot (2006), were very similar and provided the best model for the observed glint.The conclusions remain valid for (A)ATSR as it has similar channels to MODIS.Therefore the anisotropic model of Cox and Munk (1954a) is adopted here; the resulting expression for the slope distribution is where the terms ζ = Z x /σ x and η = Z y /σ y , and σ x and σ y are the root mean square values of Z x and Z y , respectively.The values σ 2 x , taken as as 0.003+0.00192w±0.002, and σ 2 y , taken as 0.00316 w±0.004, are from Cox and Munk (1954a) for a clean sea surface.
The total contribution ρ gl is calculated, following Cox and Munk (1954a), as www.atmos-meas-tech.net/3/813/2010/Atmos.Meas.Tech., 3, 813-838, 2010 where R f is the Fresnel reflection coefficient, and β the facet tilt where the scattering angle between the incident beam and surface facet is obtained using such that 2 is the familiar scattering angle between incident and reflected beams.Calculation of the Fresnel reflection coefficient R f requires the real component of the refractive index of air, taken as n a =1.00029 for all wavelengths, and of water n w , shown in Table 3.These were calculated at 550 nm and 670 nm using the method of Quan and Fry (1995) assuming a temperature of 15 • C and salinity of 35 parts per thousand, but are correct to four significant figures over the range of typical temperatures and salinities.
This model of Quan and Fry (1995) extends only to 700 nm, so for the longer wavelengths values for pure water from Hale and Querry (1973) were used.At shorter wavelengths there was an offset of around 0.0065 between the refractive index as predicted for pure water and that of salinity typical for the sea, so this adjustment was also applied to the pure water data used at 870 nm and 1.6 µm.Values for the imaginary part of the refractive index were likewise taken from Hale and Querry (1973).
For viewing zenith angles more extreme than 70 • , a slight modification is made to the denominator of Eq. ( 12) following Zeisse (1995) to avoid an infinite radiance at viewing zenith angles near the ocean horizon.The reader is referred to Zeisse (1995) for more details; although such extreme viewing geometries do not occur for the (A)ATSR views, calculation of the reflectance at such geometries is required for the integration to obtain R bd and R dd for the retrieval forward model.

Magnitude of contribution
The glint contribution ρ gl to R bb is shown at 550 nm in Fig. 4. The strong geometric dependence is visible; this is key to the ability of (A)ATSR to perform retrievals into the sunglint region, as generally while one view is affected by glint (meaning most signal arises from the surface) the other is not (so most signal arises from the atmosphere).The asymmetry in Fig. 4 arises due to the wind direction not being in line with the field of view.This has a smaller impact as θ v tends to the nadir.Dependence on wavelength is weak due to the similarity of the refractive index of water at the modelled wavelengths.
The sea surface BRDF is integrated using Gaussian quadrature with 4 points to obtain R bd and R dd .The glint contribution requires a large number of points to calculate accurately; as a result precalculated lookup tables (LUTs) of integrated ρ gl , using 360 quadrature points, are used for computational efficiency.These are parametrised in terms of wind speed (for R dd ) and wind speed and solar zenith angle (for R bd ).These quantities are shown in Figs. 5 and 6.The glint DHR shows the expected increase with solar zenith angle, and for commonly-encountered conditions is approximately 0.03 (at all wavelengths; shown only for 550 nm).For high wind speeds the contribution decreases due to the increased whitecap fraction.The BHR is generally 0.05-0.06,decreasing as w increases again due to the increase in whitecap fraction.When the whitecap contribution is added, an increase with w is observed (except at 1.6 µm where the foam reflectance is small) together with more variability between the wavelengths.

Uncertainties
Over a range of typical conditions, the uncertainties in the coefficients used in the calculation of p(Z x ,Z y ) (Eq. 11) lead to a variability of around 10%, causing a corresponding uncertainty in ρ gl of the same amount.This variability decreases at higher wind speeds.

Underlight
Underlight is upwelling irradiance from just below the surface of the ocean.As such, ρ ul is influenced strongly by pigment concentration and wavelength, and weakly by geometry.The model described here is designed for Case I waters, following the nomenclature of Morel and Prieur (1977).In Case I waters (typically open ocean) the chlorophyll concentration is high compared to the scattering coefficient; in Case II waters (typically coastal and shallow) scattering by inorganic particles dominates.The semi-empirical relationships between ocean constituents and surface reflectance es.developed by Morel and Prieur (1977), Morel and Gentili (1991) and later work are different between Case I and II waters.This work focusses on Case I waters because they cover the majority of the Earth's surface, and scattering in Case II waters is less well-understood.

Calculation
The underlight reflectance is an analogue to the atmospheric scattering problem.Reflectances and transmittances related to underlight are denoted using R and T , rather than R and T , for an easier distinction between other reflectance and transmittance terms used in this work.Fundamentally, the system may be considered to consist of three layers:   5), and the solid lines the total contribution from whitecaps and glint.Fig. 6.Glint and whitecap contribution to bihemispherical reflectance R dd as a function of wind speed.The black lines shows the contribution at 550 nm, red at 660 nm, green at 870 nm, and the blue line the contribution at 1.6 µm.The dashed lines show the contribution from sun-glint alone, weighted against by the whitecap fraction (Eq.5), and the solid lines the total contribution from whitecaps and glint.
-The topmost layer, corresponding to the air-water interface.The downward and upward transmittances through this surface are denoted T d and T u , respectively; the (downward) reflectance of upwelling irradiance below the interface is R u .The upward reflectance of downwelling irradiance above the interface is the glint.
-An "ocean floor" layer.This is characterised by the reflectance of the ocean bed, R bed .The transmittance between the upper layer and ocean floor is denoted T w,d or T w,u for downward and upward directions, respectively.
The problem may be simplified by first considering a combination of the lower two layers.Light penetrating the airwater interface may be reflected back towards it (R w ) or subject to multiple "reflections" between the upper ocean and ocean floor.If it is assumed that R bed is isotropic then it can be easily shown that the combination of these two layers reduces to which is the reflectance R w of the incident light from the water body, plus a multiple-scattering geometric series limit.This combination of the two lower layers may then be treated as a single (lower) layer of a two-layer system, in which the upper layer corresponds to the air-water interface.If it is assumed that this lower layer is an isotropic reflector then the same series limit may be applied to this simpler two-layer system to calculate ρ ul : The geometric dependence of R w is weak and its absolute value is small, so the error introduced by this approximation is small.The direct reflectance of incoming light off the interface was dealt with as the sun-glint term and so is not part of Eq. ( 16).A further approximation may be made to simplify the calculation of ρ ul for the open ocean.The transmittance T w of water (either upward or downward) may be calculated as where a w is the absorption coefficient of the water, and z the path length (for a vertical column, equal to the depth of the water).For pure water, a w can be calculated from the complex part of the refractive index: Using κ from Table 3, over all wavelengths of interest, and even in shallow water (with z = 100 m), T w is very small (10 −2 for z = 100 m at 550 nm, and orders of magnitude smaller for deeper water or longer wavelengths).Dissolved substances in seawater would further decrease T w .As a result T w,d R bed T w,u , the proportion of light transmitted through the water, reflected off the bottom and then transmitted up through the water body, is negligible and Eq. ( 16) may be simplified to which is the expression used to calculate ρ ul in this scheme.
It is noteworthy that in very shallow waters, or wavelengths at which water is more transparent, the reflectance characteristics of the ocean floor may become important and so Eq. ( 16) is presented as the general case.An analagous formulation was presented by Austin (1974).

Downwelling transmittance coefficient, T d
The term T d in Eq. ( 19) represents the transmittance of downwelling radiation.Assuming a flat sea surface, this is simply calculated using the Fresnel coefficient for an incident beam of a given solar zenith angle, and noting that light not reflected is transmitted: The subscript in R f:aw reminds that the incident beam is coming from the air, into the water.For all wavelengths, T d is approximately 0.98 for θ s <60 • but drops sharply for larger zenith angles.Calculation for a wind-roughened sea is computationally expensive, as it involves the calculation of the transmittance through all possible facets.Austin (1974) present results for selected angles and wind speeds, and note that for wind-ruffled seas T d is slightly lower than 0.98 for near-zenith angles of incidence, and the decline in transmittance is slower as the Sun approaches the horizon, although the changes are not large.Therefore the assumption of a flat sea surface introduces minimal additional error.

Upwelling transmittance coefficient, T u
The transmittance of the underlight through the water-air interface is denoted T u .If the upwelling irradiance is assumed to be diffuse, and the sea surface flat, then T u is given by integrating the Fresnel coefficient over all possible upwelling angles θ u : Here, R f:wa indicates that the upwelling light is travelling from water to air.This gives T u =0.522 at 550 nm, 0.523 at 660 nm, 0.525 at 870 nm and 0.536 at 1.6 µm.These are just over half typical values of T d because rays hitting the interface with > sin −1 (n a /n w ), approximately 48 • , are internally reflected so that their energy is lost.Hence the radiance penetrating the surface is limited to the subset with angles of incidence smaller than this critical angle.
As with T d , for a rough sea calculation becomes more complicated because the transmittance of facets aligned at different angles to the surface has to be taken into account.Austin (1974) again show, using examples at selected wind speeds and angles, that for increasingly rough seas the transmittance from upwelling rays at near-nadir incidence falls while some transmittance is possible for rays at angles larger than the flat-sea critical angle.The net effect is that T u shows little dependence on wind speed, and so the flat-sea assumption is again valid.

Upwelling reflectance coefficient, R u
The final geometric term R u is the (downward) reflectance coefficient for upwelling radiance at the water-air boundary.This can be calculated as 1 − T u .Austin (1974) give broadband visible values between 0.485 for a still ocean surface and 0.463 for a wind-ruffled surface with w = 16 ms −1 ; as this dependence on wind speed is small, and R u R w 1, the flat-surface assumption introduces negligible error into Eq.( 19).

Water body reflectance, R w
The water body reflectance R w is controlled by the optical properties of water and matter within it, and is defined as the ratio of upwelling irradiance from just below the surface to downwelling irradiance just above it.The method of calculation is based on the method of in Morel and Prieur (1977), and further developed on many occasions (e.g. in Morel, 1988 or Morel andGentili, 1991).The parametrisations are based on a variety of semiempirical relationships.The water body reflectance is calculated from the optical properties of the water as follows: This describes the colour of the water as the ratio of the total backscattering coefficient b b (λ) to the absorption coefficient a(λ), multiplied by some empirical correction factor f .

Absorption coefficient
A more thorough treatment can be given to the absorption coefficient of water than the approximation made previously.The total absorption coefficient a of seawater can be thought of as the sum of the absorption due to pure water, a w (as in Eq. 18), that due to phytoplankton pigments a ph , and a CDOM , the absorption due to detritus and coloured dissolved organic matter (CDOM), also known as Gelbstoff : The absorption coefficients used for water are shown in Table 4. Values for 550 nm and 660 nm are taken for seawater from Morel and Prieur (1977); for longer wavelengths data are unavailable so at 870 nm and 1.6 µm a w is estimated using imaginary components of the refractive index from Table 3 with Eq. ( 18).This approximation is justified as the underlight contribution to R bb is small at these wavelengths.
For a ph , the two-component model outlined by Sathyendranath et al. (2001) and Devred et al. (2006) is used.This relates the absorption due to phytoplankton to the concentration C of chlorophyll-a in mg m −3 , assuming a mixed population of two phytoplankton types, by where the parameter U is defined as where C m 1 is the maximum chlorophyll-a concentration associated with phytoplankton population 1 in mg m −3 , and a * 1 and a * 2 are the specific absorption coefficients in m −1 (mg chl-a) −1 of the two populations at the wavelength of interest.S chl describes the nonlinearity of absorption and has units of m 3 (mg chl-a) −1 .For a global dataset, Devred et al. (2006) found C m 1 = 0.62 mg m −3 , a * 1 = 0.0109 m −1 (mg chla) −1 at 550 nm and 0.0173 m −1 (mg chl-a) −1 at 660 nm, a * 2 = 0.0064 m −1 (mg chl-a) −1 at 550 nm and 0.0085 m −1 (mg chl-a) −1 at 660 nm, and S chl = 1.61 m 3 (mg chl-a) −1 .Absorption by pigments is neglected at 870 nm and 1.6 µm; the very strong absorption of the water at these wavelengths (Table 4) means this approximation has negligible impact.
For pigment concentrations of approximately 1 mg m −3 or more a ph becomes significant at 550 nm and 660 nm, otherwise it is negligible compared to a w .Operationally, data for both chlorophyll concentration and CDOM/detritus absorption are obtained from the Glob-Colour project (Barrot et al., 2006).This provides global values of various ocean colour parameters from merged satellite (MERIS, SeaWIFS and MODIS) datasets.Monthly mean values on an approximately 25 km×25 km grid are used, with gaps filled using an annual mean.Figure 7 shows the annual mean pigment (chlorophyll-a and phaeophytin-a, normally abbreviated as just "chlorophyll") concentrations over 2004.The large spatial variability, as well as range of concentrations spanning several orders of magnitude, is evident.In the open ocean, typical values are in the range 0.05 to 1 mg m −3 but near the coast the concentration can reach 10 mg m −3 or higher.Chlorophyll content is also important for calculating the backscattering coefficient, as discussed in the next section.
Figure 8 presents an analagous map of the annual mean CDOM absorption coefficient at 550 nm for 2004.According to Roesler et al. (1989), absorption from detritus and CDOM can be treated as one parameter due to their similar spatial distributions and absorption properties; hence, the quantity retrieved by ocean colour algorithms is the total absorption coefficient for both substances.GlobColour provides the absorption coefficient at 443 nm which is related to absorption at longer wavelengths by the following equation: In Eq. ( 26) the parameter S describes the spectral slope of the absorption.Roesler et al. (1989) found different values (generally in the range 0.011 to 0.018) worked well for different regions of the world, with 0.014 a good value for global studies.This value of 0.014 is used here, as well as assorted other studies (such as Chen et al., 2003).
The CDOM absorption coefficient tends to covary with chlorophyll concentration (Figs. 7 and 8).Typical ocean values at 550 nm are in the range 0.001 to 0.01 m −1 but again higher values, generally up to 0.1 m −1 , can be observed in productive or coastal waters.This algorithm only takes a CDOM into account at 550 nm.At 660 nm the value of S means that CDOM absorption is only around a fifth as strong as at 550 nm.Combined with the fact that absorption by water and chlorophyll increases by roughly an order of magnitude, the CDOM contribution to the total absorption coefficient is negligible.At longer wavelengths this effect is even more pronounced.Taking into account the decreasing importance of ρ ul with increasing wavelength, this approximation has minimal effect on results.

Backscattering coefficient
The term b b (λ) is the total backscattering coefficient from molecules and particles.Following Morel and Prieur (1977) this is parameterised as where b w is the molecular scattering coefficient for water, and the second term is the product of the particle backscattering probability bb and particle backscattering coefficient b.The division by 2 of b w arises because molecular scattering is forward-back symmetric, so the backscatter coefficient is half of the molecular scattering coefficient.
Values for b w for pure water from 380 nm to 700 nm were given by Morel and Prieur (1977).Morel (1974) tabulated values from 350 nm to 600 nm for both pure water and typical seawater.The data were shown to fit a power law with a dependence on λ −4.32 , with seawater scattering around 1.30 times as much as pure water.This relationship has been used to extrapolate these data to the 660 nm, 870 nm and 1.6 µm channels.The values obtained for b w are 1.93×10 −3 m −1 at 550 nm, 8.77×10 −4 m −1 at 660 nm, 2.66×10 −4 m −1 at 870 nm and 1.91×10 −5 m −1 at 1.6 µm.The value predicted for 660 nm is in good agreement with that given for pure water at 660 nm in Morel and Prieur (1977) multiplied by 1.30.The small size of b w (both in absolute terms and when compared to bb ) at longer wavelengths means that any error in this extrapolation is minor in terms of influence on R w .
The second parameter in Eq. ( 27), bb (λ), is the backscattering probability: the ratio of the backscattering to scattering coefficients of the pigments.It is related to the total concentration C of chlorophyll=a and pheophytin-a, measured in mg m −3 , and wavelength λ, measured in nm, by the following expression: bb (λ) = 0.002 + 0.02(0.5 − 0.25log 10 C) 550 λ (28) The final term in the backscatter component of Eq. ( 27), b is calculated as: b = 0.3C 0.62  (29)   The relationship between b and C was derived by Morel (1988) for data at 550 nm; the wavelength-dependence of particle backscattering is taken into account by the λ −1 factor in Eq. ( 28).It should be noted that although parametrised in terms of C, the models were developed to account for scattering from suspended organic matter as well as pigment.

Ratio multiplier f and combination for water
body reflectance R w Morel and Prieur (1977) initially gave f , the empirical correction multiplier of the ratio of total backscattering to total absorption used to calculate the water body reflectance R w , a value of 0.33.Subsequent work has found it to depend on the solar geometry and the optical properties of water.The method used here was put forward by Morel and Gentili (1991), stated to be accurate within 1.5% for solar zenith angles smaller than 70 Assuming a pigment concentration of 0.3 mg m −3 , the variation of f is small with wavelength but larger with solar angle, from slightly over 0.3 for a near-nadir sun to over 0.5 for a sun near the horizon.

Magnitude of contribution
Figure 9 shows ρ ul at (A)ATSR wavelengths for a range of representative pigment concentrations.At the shorter wavelengths it is of the order of 10 −2 − 10 −3 , and so away from the glint region is generally equal to or larger than other contributions to R bb .Hence knowledge of C is essential to judge accurately the total reflectance.As it shows a stronger wavelength-dependence than ρ wc and ρ gl , the spectral shape of R bb will be largely determined by ρ ul outside of the sunglint region.At 870 nm and 1.6 µm, ρ ul is negligible.At all wavelengths ρ ul increases with C, although at 550 nm ρ ul decreases for C > 1 mg m −3 as the increased a CDOM(550) used in the calculations causes absorption to increase more rapidly than scattering.

Uncertainties
The major uncertainties associated with ρ ul are errors arising from poor characterisation of pigment and CDOM distributions and scattering.The small size of the underlight term at the longer two wavelengths means that errors in ρ ul will only have minor impacts on the modelled reflectance at 550 nm and 660 nm.
-The relationship between bb and C was developed for Case I waters (according to the definitions of Morel and Prieur, 1977) and so may not accurately characterise scattering in Case II waters (where pigment and scattering particles do not covary in the same way).This may cause the algorithm to perform less well over Case II waters.Case II waters are largely coastal and the inhomogeneity of coastal regions presents other problems for aerosol retrieval; this is beyond the scope of this work.
-Errors arising from use of monthly means for chlorophyll and CDOM values.The GlobColour chlorophyll products have a stated accuracy of 31%.CDOM errors are not given by Barrot et al. (2006).Further errors arise due to variations on shorter timescales than a month.
-There may be regional biases from using 0.014 as a global CDOM spectral slope S, as Roesler et al. (1989) found values from 0.011 to 0.018 in different parts of the world.

Typical values, patterns, and errors
To illustrate typical reflectances predicted by the model, it has been applied to a selection of AATSR orbits from 6 September 2004.Information on the viewing solar-satellite geometry is shown in Fig. 10, along with the ECMWF nearsurface wind speed from noon on that date.Figures 11 and 12 show example nadir-view and forwardview BRDFs generated at 550 nm and 1600 nm.Sun-glint is visible in the centre of the nadir-view swaths and toward the north of the forward-view swaths.These glint patterns persist for different orbits as the satellite geometry remains the same.Perturbations to this glint shape arise due to variations in the wind speed and from increased reflectance due to oceanic whitecaps or underlight, particularly at 550 nm.BRDF at 660 nm and 870 nm takes values in between these two wavelengths; at 550 nm all contributions to the BRDF are important, while at 1600 nm the shape is glint-dominated, unless wind speeds are high.The exact location of the glint region varies seasonally.
The sea surface BRDF for (A)ATSR is observed to take values from around 10 −5 to 1 dependent on the wavelength and location with respect to the sun-glint region.The highest values are observed for shortest wavelengths; near the glint region the BRDF can be nearly spectrally flat while far away there can be orders of magnitude difference.
The DHR, shown at 550 nm and 1600 nm in Fig. 13 for the same swaths as the previous BRDF example, is almost identical for both of AATSR's viewing geometries (whose solar angles differ by under 1 • for any given pixel).Values increase as the solar zenith angle increases; as with the shape of the BRDF, perturbations to the basic shape arise due to the wind and pigment distribution.The same scale is used as in FigS.11 and 12 to illustrate the comparative variability of the reflectances.Again, values at 660 nm and 870 nm are intermediate between these two wavelengths.The DHR typically takes value around 10 −2 (where the Sun is high) to 10 −1 (where the Sun is low), and is slightly larger at shorter wavelengths.
Figure 14 shows the a priori BHR generated for these overpasses.Being independent of geometry it is the same for both instrument viewings and variability over the globe is determined by wind and pigment distributions: the glint shape of the BRDF is "averaged out" by the integration, which contributes around 0.05 to the albedo.The BHR is spectrally flatter and less variable than the BRDF or DHR.Typical values are in the range 0.05-0.08 at all wavelengths, with the shortest wavelengths being brightest.The same scale is again used as in Figs.

A priori uncertainty
It is important to assign a reasonable error to the a priori albedo generated: too small a variability and the retrieval will be unduly constrained by an imperfect model, but too large a variability and some of the information on the state is effectively thrown away.Appropriate uncertainties have been determined in the following way: -Generation of 10 000 random sets of typical ocean and viewing states (for example differing wind speeds, chlorophyll concentrations and geometries).
-For each ocean state, generation of R bb and integration for R bd and R dd for an ensemble (50 members each) of random perturbations to the uncertain model parameters (such as the foam reflectance, or sea slope characteristics from Cox and Munk, 1954a).The magnitude of the perturbations is determined by the stated uncertainty on the model parameter as previously described in the text.
-Calculation of the ensemble median reflectance and its standard deviation for each ocean state, for each of R bb , R bd and R dd .Use of medians decreases the sensitivity to outliers.
Figure 15 shows the calculated median and standard deviation of R dd for the ensembles considered.The standard deviation is observed to vary with the magnitude of R dd .The ratio of the ensemble standard deviation to the median is a measure of the proportional sensitivity of the state to errors in  the model parameters: a larger ratio means that the calculated reflectance or albedo is more sensitive to errors in the model parameters.Over all states, the median of these sensitivity ratios for R dd is 0.20 to 2 decimal places at all wavelengths.The similarity between wavelengths is an indicator that the dominating terms (the integrated glint contribution), are similar for each.In the retrieval, multiplying these sensitivities by the a priori BHR gives the a priori error estimate for the BHR for each ground scene.
As the ratios R dd : R bb and R dd : R bd are fixed in the retrieval forward model, analagous values for the sensitivity of R  ratios.For the DHR, these median relative uncertainties are 0.22 at 550 nm and 0.20 at longer wavelengths.As the BRDF spans several orders of magnitude, dependent on wavelength and position relative to the sun-glint, calculation of a relative uncertainty is not always useful (as the relative uncertainty on a dull BRDF may be high, even when the abso-lute uncertainty on it is not).As a result some minimum threshold for the uncertainty is imposed based on the median absolute uncertainty determined for dull BRDFs.The median relative uncertainty on R bb for the nadir view is 0.81 at 550 nm, 0.75 at 660 nm, 0.69 at 870 nm and 0.63 at 1.6 µm with minimum absolute uncertainties of 0.01, 0.008, 0.006   and 0.005, respectively.For the forward view, relative uncertainties are 0.82, 0.73, 0.64 and 0.58 with minimum values of 0.007, 0.004, 0.002 and 0.001 for the four channels.The higher uncertainties at shorter wavelengths arise due to the fact that whitecaps, glint and underlight may all contribute significantly to R bb at these wavelengths.
For more details on the calculation of the contribution of uncertainties in the fixed ratios R dd : R bb and R dd : R bd to the forward model error budget the reader is referred to Sayer (2008).The calculated contribution in terms of percentage uncertainty on the TOA reflectance is given in Table 5.These values are of a similar order of magnitude to the measurement uncertainty (Smith et al., 2001(Smith et al., , 2008)), and increase with wavelength as the atmospheric contribution to TOA reflectance decreases.Forward-view values are lower due to the longer atmospheric path length.The corollary of this is that, when one view is observing the glint region (bright TOA reflectance), its measurements receive less weight in the retrieval.

Application to aerosol retrievals
The new sea surface reflectance algorithm has been used to peform aerosol retrievals from AATSR data in the southeast Pacific (60 • S-0 • S; 180 • W-60 • W) for the month of September 2004.This region is chosen because it contains a large region of open ocean, far from strong aerosol source regions.The background mid-visible AOD in these open oceanic regions is typically <0.1 (Smirnov et al., 2009), meaning an accurate model of the surface reflectance is of importance to determine the surface and atmospheric contributions to the TOA reflectance.Additionally, this area contains coastal regions of high chlorophyll concentration and open regions of low chlorophyll concentration (Fig. 7), and while wind speeds over the bulk of the area are from 5-8 ms −1 , the region from approximately 45 • S-60 • S typically has stronger winds in the region of 10 ms −1 or higher (see Fig. 3).For these retrievals, the most recent version of the AATSR visible channel calibration trend data are used (v12).Results are only presented for sea retrievals.
The OE method (Rodgers, 2000) allows an analysis of retrieval performance through examination of retrieval statistics.Additionally, comparisons of aerosol optical depth with satellite data from the MODIS-Terra instrument (Remer et al., 2005), and with ground-based observations from MAN cruises in the region (Smirnov et al., 2009), and the AERONET site on the island of Tahiti (Holben et al., 1998), are possible.

Aerosol retrieval statistics
The principle of OE is to maximise the conditional probability of the retrieved state given measurements and any a priori information.Formally, this is the maximum of P = P (x|y,x a ,b) with respect to the values of the state vector x for a measurement vector of reflectances y, where x a is the a priori value of the state vector and b are all other parameters not modelled by the forward model.The assumption is made that errors in the measurements, a priori and model parameters have Gaussian distributions with zero mean and covariance matrices given by S y , S x and S b , respectively.Following Rodgers (2000) the maximum probability is given for the minimum of J , the retrieval cost: The terms present in the equation represent weighted deviations from measurements and the a priori state.Here y(x) refers to the values of y predicted by the forward model from the current value of the state vector; for clarity, the measurement vector is denoted by y m .The minimisation is done with respect to x, so that the derivative of J is independent of b.The impact of S b , the model parameter error, on J is included by mapping it into measurement space and including it as a contribution to S y .In this case examples of these model parameters include the fixed ratios R dd : R bb and R dd : R bd , set by the surface reflectance model.
Operationally, J is normalised by the number of measurements (here, 8) before being output by ORAC.Standard quality controls are applied to retrievals to exclude those poorly fit (typically a result of cloud contamination).This involves considering only those scenes retrieved with a normalised cost lower than 10.Additionally, only the maritime aerosol model was used.The distribution of residual costs is shown in Fig. 16, for both the total cost and the components corresponding to deviations from measurements and a priori.Retrieved AOD and effective radius have a large a priori uncertainty, so J is dominated by contributions from the fits of the measurements, and retrieved minus a priori BHR.
With 8 measurements and 4 constrained state vector elements, it is expected that a well-fit retrieval will have normalised cost on the order of 1.5 (=12/8), with approximately 1 coming from measurements and 0.5 from a priori.The overall distributions should correspond to χ 2 distributions with these numbers of degrees of freedom.There should be few cases of costs exceeding triple these values (approximately 4.5 for the total cost).Figure 16 shows that the distribution of actual residuals broadly matches these theoretical considerations, suggesting the uncertainties in the retrieval are well-characterised.The total cost is a good match for a χ 2 distribution for 1.5 degrees of freedom.The figure is truncated at J = 5; the number of retrievals with J > 5 is negligible.Slightly higher proportions of retrievals with 0.5 < J < 1 than expected are likely due to the fact that there are constraints, albeit weak, from the a priori AOD and effective radius so the true number of degrees of freedom is actually slightly higher than 1.5.Due to the low number of retrievals with 4.5 < J < 10, the results presented here do not change significantly if a stricter cost threshold is adopted.
As discussed, when the residuals (y(x) − y m or x − x a ) are weighted by their uncertainties (square roots of appropriate elements of S y or S x ) these distributions should each approximate Gaussians with mean 0 and variance 1.As well as the total cost, the distributions of error-normalised residuals may be examined individually to see whether each measurement or parameter is particularly well or poorly fit. Figure 17 shows the normalised residuals on the measurements, and Fig. 18 on a priori BHR.It would be expected that approximately 68% of the data should fall within the range ±1 and 99% within ±3.These figures show that, on the whole, the distributions meet these expectations, which shows again a good representation of uncertainties in the retrieval algorithm.Some distributions show a degree of bimodality, which may indicate differing deficiencies in the aerosol or surface models in some situations, although biases are small.The residual distributions on white-sky albedo are slightly narrower than expected, indicating that the BHR predicted by the reflectance model is slightly more precise than the error analysis estimated.The 550 nm residual distribution has a small positive bias, suggesting that surface reflectance at this wavelength is slightly brighter than the model predicts.
The exception is the 1.6 µm reflectance residual, which tends to be negative and wider than expected for the nadir view.As AOD is generally much lower in the nIR than visible region, the bulk of the information on optical depth is obtained from the shorter-wavelength channels, and so the poorer nadir-view 1.6 µm fit should not lead to significant errors in AOD.The width suggests that, in some situations, the forward model error for this view and channel should be larger to represent accurately the precision with which this measurement can be fit.More generally, the measurement residual distributions have a larger negative tail for the nadir view and positive for the forward view.The reasons for this are uncertain and may be due to either a poorer representation of aerosol or surface properties at nadir-viewing geometries, or alternatively issues with the data put into the surface reflectance model.An additional possibility is calibration issues over dark targets; existing vicarious calibration of the ATSR instruments has focussed on bright targets (Smith et al., 2002(Smith et al., , 2008)).A small calibration offset may be hidden in the signal from a bright target, but count more for an otherwise dark signal.

Intercomparison of aerosol optical depth
Retrieved AOD from ORAC-AATSR is compared with that obtained from the MODIS sensor aboard the satellite Terra.ENVISAT's daytime overpass is approximately 10:00 a.m.local solar time; the Terra platform shares a similar overpass time of 10:30 a.m., as opposed to Aqua's 01:30 p.m. From this point, MODIS will be taken to refer to MODIS-Terra.The QA-weighted mean ocean AOD and standard deviation, from the Collection 5 level 3 daily MODIS atmosphere product (MOD08 D3), are used.This is provided on a 1 • grid; before the comparison, the AATSR retrievals are aggregated onto the same grid.AATSR retrievals are weighted by the relative uncertainty on retrieved AOD, as provided by OE.At least 6 retrievals must be present in each bin.To migitate the effects of sampling (MODIS's swath is over 2000 km as compared to the approximately 550 km AATSR swath, and they are not on the same orbital track), the gridded data are compared only when both MODIS and AATSR provide aerosol retrievals on a given day.
Monthly mean fields of 550 nm AOD have been calculated from the daily mean fields for both instruments, and are shown in Fig. 19.Coverage is incomplete and largely limited by AATSR's narrower swath, MODIS elimination of high-sediment and sun-glint regions and overall high cloudiness.This low sampling means that features in individual orbits may still be seen in the monthly means.There is a good spatial correlation between the two, although AATSR retrieves lower AOD (typically 0.02-0.1)than MODIS (typically 0.07-0.2).Elevated AOD in open ocean regions in both correspond to higher aerosol loadings from windier conditions; although various relationships between wind speed and marine AOD have been proposed (for example, recently Smirnov et al., 2003, Glantz et al., 2009;Huang et al., 2009;Lehahn et al., 2010), all suggest higher winds lead to higher AOD.
Aerosol data from cruises in the Maritime Aerosol Network are available for the study region.The low spatial and temporal coverage of these limits the ability to directly compare with satellite data, although they may be used to provide some information about typical AOD for different regions and its variability.The cruises used are as follows: -NOAA Ronald H. Brown Cruise, December 2007-Feburary 2008.Latitudinal transect through the eastern Pacific.
-NOAA Ronald H. Brown Cruise, October-November 2008.Measurements from the vicinity of Arica, Bolivia out to the west, then north towards the Ecuadoran coast.
Across the Pacific coast of South America.
-RV Hesperides Cruise, February-March 2009.Across the Pacific coast of the southern part of South America.
Transect from Brisbane, Australia to Valparaiso, Chile.
Measurements during these cruises were made using Microtops sun photometers (Smirnov et al., 2009); although AOD at 550 nm is not retrieved from these measurements, it is estimated from the AODs recorded at 440 nm and 870 nm and the Ångström exponent between these two wavelengths.Level 2.0 data (cloud-cleared and quality-assured) were used for all cruises except RV Melville 2009-2010, where level 2.0 was not yet available so level 1.5 (cloud-cleared and corrected for pointing errors) were used.The daily average AOD at 550 nm reported from these cruises is shown in Fig. 20.The standard deviation about these daily values was generally small (<0.01).Highest optical depths, on occasion exceeding 0.4, were found in coastal regions.Typical AOD for the open ocean is in the region 0.05-0.09,interim between χ 2 disribution for total cost.MODIS and AATSR.Where both are available, in most cases the difference in daily mean AOD between the level 1.5 and level 2.0 data is negligible.For the few cases in which the level 1.5 and 2.0 data are significantly different, level 1.5 data generally have AOD>0.2,and the corresponding level 2.0 AODs are lower by 0.01-0.04.

Nadir-view fit residuals
The AERONET site at Tahiti (17.577 • S, 149.06 • W, 90 m elevation) did provide aerosol measurements during September 2004; these are available at level 1.5, and here 550 nm AOD is estimated from that retrieved at 500 nm, 870 nm and the Ångström exponent between these wavelengths.A time series of daily mean 550 nm AOD is shown in Fig. 21, with the standard deviation about the daily average given as an uncertainty estimate.Also shown are the MODIS and AATSR mean and standard deviation AOD for the 1 • grid box in which Tahiti falls.The 550 nm AOD measured at the AERONET site is typically between 0.03-0.07.Although the number of cloud-free satellite overpasses of the region is low, the figure suggests a tendency for MODIS to overestimate the AOD as compared to AERONET (although the variability of the MODIS data are high), and AATSR to underestimate.Disagreement would also arise if the Tahiti AERONET site were not representative of the larger area, although in this case the impact would be the same for both satellite datasets.
Histograms of the monthly mean satellite AOD, as well as the ground-based measurements, are presented in Fig. 22.As the majority of the satellite observations are open ocean (Fig. 19), the MAN cruises are restricted to those two which went through this region, Ronald H. Brown (2007Brown ( -2008) ) and RV Melville (2009Melville ( -2010)).Additionally, the AOD observed at the AERONET Tahiti site is shown.All ground-based measurements, as opposed to daily averages, were used to generate the histograms.Together, these give some indication of typical AODs found in the region, although neither satellite dataset would be expected to match as these may not represent the whole region well.Also shown is a histogram corresponding to the AATSR data reprocessed assuming a fixed wind speed of 6 ms −1 , and the overall surface reflectance fixed at the value predicted by the model, as in the MODIS Collection 5 algorithm (Remer et al., 2005).cessed, the overall distribution of AOD is observed to slightly broaden and the peak shifted to higher AOD by around 0.01.This indicates that some of the broadness and positive offset of MODIS as compared to AATSR may be explained by the assumption about a fixed surface reflectance and wind speed.In this respect, the assumptions about wind speed made in the MODIS aerosol retrieval are likely to lead to an overestimate of AOD in windy (w > 6 ms −1 ) conditions.A higher MODIS AOD could also arise if MODIS underestimated the underlight reflectance, although this is small over large parts of the region considered due to low chlorophyll a concentrations, and so unlikely to be important.Kaufman et al. (2005) examined the effect of residual cloud contamination on MODIS AODs over ocean and concluded that at 550 nm contamination, principally from thin cirrus, could  lead to a false enhancement in AOD of the order of 0.02 at 550 nm.An analagous analysis has not been carried out for AATSR, although the cost function has been generally found to be successful at identifying cloud-contaminated scenes.Together these explain some of the discrepancy between MODIS and AATSR AODs.Kokhanovsky et al. (2009) present an intercomparison study of aerosol retrievals over a black surface performed from synthetic data.Both the ORAC-AATSR and MODIS ocean aerosol retrieval algorithms participated in this study.In Kokhanovsky et al. (2009) was observed to have a positive bias and ORAC-AATSR a slight negative bias in AOD.At low optical depths, as observed in this work, these biases were small (of the order 0.01 at 550 nm).Additionally, it is not clear whether the retrievals using idealised synthetic data reflect accurately the performance of the retrieval using real measurements.
Calibration differences between the two instruments over these dark targets, and assumptions made about aerosol models, may also lead to disagreement.MODIS observes larger AODs on average than either the Tahiti AERONET site, or the MAN cruises, although this may reflect the conditions at these sampled locations, or biases in the sun-photometer retrieval algorithms.The same factors may explain the lower AODs as seen by AATSR as compared to the ground-based measurements.It is also possible that AATSR is overestimating the surface reflectance, and thereby underestimating the AOD.However, as shown in the previous section, the retrieval costs indicate that the retrieved states are consistent with the measurement and a priori uncertainty, with residuals reasonably unbiased (except in the nadir view at 1.6 µm), so it is difficult to diagnose a possible bias in retrieved AOD by examining fit statistics.

Conclusions
A sea surface BRDF model, drawing on the heritage of Koepke (1984), has been formulated for the visible and near-IR channels of the ATSR instruments (ranging from 550 nm to 1.6 µm).The model accounts for contributions to the observed reflectance from whitecaps, sun-glint and underlight.The model is discussed in the context of application to aerosol retrievals, although it is also suitable for use in cases of optically-thin clouds.It can be integrated over solar and viewing geometries to provide the DHR and BHR additionally required for the ORAC aerosol forward model.Furthermore, as ORAC is an optimal estimation algorithm, uncertainties in the parametrisations used in the BRDF model are propagated through into the forward model and retrieved state.As the brightness of the surface is permitted to vary in ORAC, unlike most other oceanic aerosol retrieval algorithms, some additional flexibility is available in the case where the assumed surface reflectance is incorrect.The new BRDF model offers improved coverage over previous methods, as retrievals are possible into the sun-glint region.
The new BRDF model has been implemented in the ORAC aerosol retrieval scheme and used to process one month of AATSR data in the south-eastern Pacific.Examination of retrieval statistics shows the assumed error budget to be generally appropriate, meaning the retrieved states are consistent with the measurements and a priori assumptions.The resulting field of AOD is intercompared with MODIS-Terra measurements in the whole region, AERONET observations at the Tahiti site, and MAN cruises over the same general area but different times.MODIS and AATSR show similar spatial distributions of AOD, although MODIS reports values which are larger and more variable.It is suggested that assumptions in the MODIS aerosol retrieval algorithm may lead to a positive bias in MODIS AOD over ocean regions where the wind speed is significantly higher than 6 ms −1 .Other differences may arise due to residual cloud contamination, calibration differences, and aerosol model assumptions.

Fig. 2 .
Fig. 2. Whitecap coverage and contribution to the BRDF as a function of wind speed.The dashed black line indicates the whitecap coverage fwc.The red line shows the contribution fwcρwc to Rbb at 550 nm and 660

Fig. 3 .
Fig. 3. Global mean near-surface wind speed from ECMWF for the year 2004, sampled at AATSR overpass times.

Fig. 5 .
Fig. 5. Glint contribution to the directional-hemispherical reflectance Rbd at 550 nm as a function of wind speed and solar zenith angle.

Fig. 5 .
Fig. 5. Glint contribution to the directional-hemispherical reflectance R bd at 550 nm as a function of wind speed and solar zenith angle.

Fig. 6 .
Fig.6.Glint and whitecap contribution to bihemispherical reflectance Rdd as a function of wind speed.The black lines shows the contribution at 550 nm, red at 660 nm, green at 870 nm, and the blue line the contribution at 1.6 µm.The dashed lines show the contribution from sun-glint alone, weighted against by the whitecap fraction (Equation5), and the solid lines the total contribution from whitecaps and glint.

Fig. 9 .Fig. 9 .
Fig. 9. Variation of ρul with chlorophyll-a concentration C at the four (A)ATSR channels used..The solar zenith angle was taken as 30 • , and a CDOM(550) (in m −1 ) was set to 10 % of C (in mg m −3 ) to represent the typically-covarying nature of these quantities in open waters as seen in Figures 7 and 8.

Fig. 10 .
Fig. 10.Viewing geometry and wind speed for AATSR swaths on 6 September 2004.The top row shows (left-right) the nadir view satellite zenith angle, forward view satellite zenith angle, and solar zenith angle.The bottom row shows (left-right) the nadir-view relative azimuth angle, the forward-view relative azimuth angle, and the 12 pm UTC ECMWF 10 m wind speed.

Fig. 10 .
Fig. 10.Viewing geometry and wind speed for AATSR swaths on 6 September 2004.The top row shows (left-right) the nadir view satellite zenith angle, forward view satellite zenith angle, and solar zenith angle.The bottom row shows (left-right) the nadir-view relative azimuth angle, the forward-view relative azimuth angle, and the 12:00 p.m UTC ECMWF 10 m wind speed.

Fig. 15 .
Fig. 15.Joint histograms of ensemble median R dd against the standard deviation of the ensemble's R dd , for (A)ATSR visible/nIR wavelengths.The colour scale indicates the relative density of points.

Fig. 15 .
Fig. 15.Joint histograms of ensemble median R dd against the standard deviation of the ensemble's R dd , for (A)ATSR visible/nIR wavelengths.The colour scale indicates the relative density of points.

Fig. 16 .FrequencyFig. 17 .Fig. 16 .
Fig. 16.Relative frequency distribution of retrieval cost, J.Total cost is shown in black, while distributions for measurement and a priori contributions are shown separately in red and green respectively.The scale indicates the proportion of all observations falling in that bin.The dashed black line indicates the theoretical χ 2 disribution for total cost.

FrequencyFig. 17 .Fig. 17 .
Fig. 17.Error-normalised residuals on nadir (left) and forward-view (right) measurements for aerosol retrievals in the southern Pacific during September 2004.In both plots, black indicates 550 nm, red 660 nm, green 870 nm and blue 1.6 µm.The vertical lines indicate 0, ±1 and ±3 respectively.The scale indicates the proportion of all observations falling in that bin.

Fig. 18 .
Fig. 18.Error-normalised residuals on a priori BHR for aerosol retrievals in the southern Pacific during September 2004.Black indicates 550 nm, red 660 nm, green 870 nm and blue 1.6 µm.The vertical lines indicate 0, ±1 and ±3, respectively.The scale indicates the proportion of all observations falling in that bin.

Fig. 19 .Fig. 20 .
Fig. 19.Monthly mean 550 nm AOD from AATSR (left) and MODIS-Terra (right) on a 1 • grid, constructed from only those grid cells where AATSR and MODIS-Terra reported observations on the same day.

Fig. 20 .
Fig. 20.Daily mean 550 nm AOD from MAN measurements in the southern Pacific.

Fig. 21 .
Fig. 21.Daily mean 550 nm AOD observed at the Tahiti AERONET site (17.577• S, 149.06 • W). Green symbols indicate AERONET measurements, and uncertainties the standard deviation of AOD throughout the day.Black triangles indicate AATSR and green diamonds MODIS-Terra, with error bars corresponding to the standard deviation of the 1 • grid cell in which Tahiti falls.

Fig. 21 .
Fig. 21.Daily mean 550 nm AOD observed at the Tahiti AERONET site (17.577• S, 149.06 • W). Green symbols indicate AERONET measurements, and uncertainties the standard deviation of AOD throughout the day.Black triangles indicate AATSR and green diamonds MODIS-Terra, with error bars corresponding to the standard deviation of the 1 • grid cell in which Tahiti falls.

Fig. 22 .Fig. 22 .
Fig. 22. Histogram of 550 nm AOD from each dataset.Black indicates AATSR and red MODIS-Terra.Both are calculated from a monthly mean on a 1 • grid, for only those cells where both instruments observe on the same day.Solid black indicates the standard AATSR retrieval, and dashed black when the wind speed is fixed at 6 ms −1 and overall surface reflectance fixed.Green indicates the Tahiti AERONET site, calculated from all observations during the month of September 2004.Blue indicates data from the open-Pacific MAN cruises Ronald H. Brown (2007-2008) and RV Melville (2009-2010), calculated from all observations taken within the study region.The scale indicates the proportion of all observations falling in that bin.

Table 1 .
Geometric notation relating to reflectance used in this work.

Table 2 .
Spectral and directional variability of components of R bb .

Table 3 .
Quan and Fry (1995)73)or water at relevant wavelengths.Both real n and imaginary κ components are given, although only the real part is used in determining the Fresnel coefficient.Adapted fromHale and Querry (1973)andQuan and Fry (1995).

Table 5 .
Forward model uncertainty in TOA reflectance, in units of the percentage of the measured signal, arising from uncertainty in the ratios R dd : R bb and R dd : R bd .