1 MODIS Ocean Science Team Algorithm Theoretical Basis Document ATBD 20 Instantaneous Photosynthetically Available Radiation and Absorbed Radiation by Phytoplankton Version 7 February 2003 Kendall L. Carder, F. Robert Chen, and Steve K. Hawes College of Marine Science University of South Florida 140 Seventh Avenue South St. Petersburg, Florida 33701 kcarder@monty.marine.usf.edu 2 TABLE OF CONTENTS 1.0 Introduction........................................................................................................................................... 3 2.0 Overview and Background Information ............................................................................................. 3 2.1 Experimental Objective.................................................................................................................... 3 2.2 Historical Perspective....................................................................................................................... 4 2.3 Instrument Characteristics................................................................................................................ 4 3.0 Algorithm Description .......................................................................................................................... 4 3.1 Theoretical Description.................................................................................................................. 4 3.1.1 Physics of Problem................................................................................................................. 4 3.1.2 Mathematical Description of Algorithm ................................................................................ 5 3.1.2.1 Calculation of Ed(8i,0+)................................................................................................. 5 3.1.2.1.1 Direct irradiance ­ Edd(8,0+) ............................................................................... 5 3.1.2.1.2 Diffuse irradiance ­ Eds(8,0+) ............................................................................. 9 3.1.2.2 Calculation of IPAR ................................................................................................... 10 3.1.2.2.1 Sea Surface Reflectance ................................................................................... 11 3.1.2.2.2 Integration of Ed over wavelength .................................................................... 12 3.1.2.3 Calculation of ARP..................................................................................................... 13 3.1.3 Sensitivity of the Algorithm ................................................................................................. 15 3.2 Practical Considerations.............................................................................................................. 17 3.2.1 Numerical Computation ....................................................................................................... 17 3.2.2 Programming/Procedural Considerations............................................................................. 18 3.2.3 Calibration and Validation ................................................................................................... 18 3.2.4 Data Dependencies ............................................................................................................... 19 3.3 Algorithm tests using SeaWiFS Data ......................................................................................... 19 3.4 Algorithm tests using MODIS Data............................................................................................ 19 4.0 Constraints, Limitations, Assumptions............................................................................................. 20 5.0 References ............................................................................................................................................ 20 6.0 Appendix -- weighting functions for IPAR and ARP calculations ................................................ 22 3 1.0 Introduction The algorithm presented here yields three related products, collectively referred to as product MOD21. The first product is the downwelling irradiance just above the sea surface in each of the visible MODIS wavebands, Ed(8i,0+), where 8i = 412, 443, 488, 531, 551, and 667 nm. This portion of the algorithm is based on the maritime irradiance model described in Gregg and Carder [1990]. The second product is instantaneous photosynthetically available radiation, IPAR, which is the total downwelling photon flux just below the sea surface, integrated over the wavelength range 400 to 700 nm. It is called "instantaneous" because it is only a measure of PAR in the instant that the sensor views a given pixel and thus does not represent the irradiance averaged over the entire day. Therefore, IPAR cannot be used directly in primary production models that require PAR values [Platt and Sathyendranath, 1988; Platt et al., 1991]. However, it may be possible to relate IPAR to daily PAR values. IPAR is most useful in measuring spatial or day-to-day differences in incident irradiance for comparison with fields of solar-stimulated fluorescence (see Dr. Mark Abbott's ATBD-MOD-23). The third and most important product is the absorbed radiation by phytoplankton, ARP. It is the total number of photons, or quanta, absorbed by phytoplankton in the top attenuation depth measured at 685 nm, z685. It is determined by multiplying the scalar irradiance and the phytoplankton absorption coefficient at each wavelength provided by MOD18 and integrating the product from 400 to 700 nm and from the surface to z685. z685 is the depth at which Ed(685,z) = Ed(685,0­)@e­1. The main use of ARP is in conjunction with the chlorophyll fluorescence algorithm (product MOD19, ATBD-MOD-23). MOD19 will provide the fluorescence line height, FLH. Dividing FLH by ARP gives a value that is proportional to the quantum yield of fluorescence, which is called chlorophyll fluorescence efficiency, CFE, in ATBDMOD-23. Even though ARP is the number of quanta absorbed by all the phytoplankton pigments, not just by chlorophyll, we will adopt the term CFE for consistency. 2.0 Overview and Background Information 2.1 Experimental Objective Each of the three products has its own experimental objective. Ed(8i,0+) is an interim product. IPAR can be used in primary production research. ARP is the most important product as it is needed to convert FLH into a value that represents the CFE of the phytoplankton. Falkowski and Kolber [1994] suggest that CFE is inversely proportional to the quantum yield of photosynthesis. Because once a photon is absorbed by a viable phytoplankton pigment, its energy must go into photosynthesis, fluorescence, or heat. While the use of FLH and CFE in estimating photosynthetic rates is the subject of much debate, the possibility of using satellites to measure primary production is enticing. CFE has also 4 been demonstrated to be related to nutrient- and/or light-limitation [Keifer, 1973a,b; Carder and Steward, 1985]. 2.2 Historical Perspective Starting with Leckner, [1978], a series of simple irradiance models have been developed, e.g., those of Justus and Paris [1985], Bird and Riordan [1986], and Green and Chai [1988]. All of these models are specific for terrigenous aerosols, which differ greatly in size and optical characteristics from marine aerosols. The total and spectral irradiance computed using these models can be quite different from the irradiance entering the ocean. The irradiance model of Gregg and Carder [1990] uses a mixture of marine and terrigenous aerosols and is well suited for maritime irradiance calculations. The Ed(8i,0+) portion of our algorithm is an adaptation of the Gregg and Carder [1990] model that uses data inputs from MODIS and other EOS sensors. Measuring global primary production is considered an important goal in oceanography. Satellite measurements of CFE may provide a means of improving estimates of global primary production (Abbott's ATBD-MOD-23). 2.3 Instrument Characteristics The bulk of the algorithm involves computations on known quantities and data products from MODIS or from other ancillary sources. The instrument characteristics important to this algorithm depend on the other algorithms. 3.0 Algorithm Description The algorithm calculates the three separate quantities sequentially, Ed(8i,0+), IPAR, then ARP. Thus, the physics and mathematics sections below will discuss each output product in turn. 3.1 Theoretical Description 3.1.1 Physics of Problem Attenuation of solar irradiance in the visible and near-UV wavelengths can be attributed to five atmospheric processes: scattering by the gas mixture (Rayleigh scattering), absorption by ozone, absorption by the gas mixture (primarily by oxygen), absorption by water vapor, and scattering and absorption by aerosols. Direct irradiance is not scattered but proceeds directly to the surface of the earth after losses by absorption. Diffuse irradiance is scattered out of the direct beam but toward the surface. The sum of the direct and diffuse components defines the downwelling surface irradiance. 5 Downward irradiance at the sea surface is then attenuated by reflection at the air-sea interface. Reflectance of the direct beam depends on the solar zenith angle and the real part of the index of refraction of seawater. Reflectance of the diffuse irradiance is related to the roughness of the sea surface. Reflectance due to foam can be related to the wind speed, and it affects both the direct and the diffuse components. The number of quanta absorbed by phytoplankton is calculated as the product of the scalar irradiance and the phytoplankton absorption coefficient integrated over the top attenuation depth. 3.1.2 Mathematical Description of Algorithm The Gregg and Carder [1990] model is an extension and simplification of the Bird and Riordan [1986] model, and the description here follows their development. The first step in the algorithm is to compute the downwelling irradiance just above the sea surface, Ed(8,0+), at 1 nm resolution. This spectrum is then binned and weighted appropriately to give the irradiance in each of the visible MODIS channels, Ed(8i,0+). Next, the below-surface values are computed, Ed(8i, Z), and summed with appropriate weights to give IPAR. Last, scalar irradiance, E0(8i, Z) is multiplied by the phytoplankton absorption coefficient, aN(8i), summed with appropriate weighting factors, and integrated over the top attenuation depth to yield ARP. 3.1.2.1 Calculation of Ed(8i,0+) Ed(8,0+) is separated into its direct and diffuse components, E d ( ,0 + ) = E dd ( ,0 + ) + E ds ( ,0 + )
where the subscripts dd and ds refer to direct and diffuse components, respectively. 3.1.2.1.1 Direct irradiance ­ Edd(8,0+) Edd(8,0+) is computed by E dd ( ,0 + ) = F0 ( ) cos( ) Tr ( ) Toz ( ) To ( ) Tw ( ) Ta ( )
where F0(8) is the mean extraterrestrial irradiance corrected for earth-sun distance and orbital eccentricity, 2 is solar zenith angle, and T is the transmittance after absorption and/or scattering by each atmospheric component. The components r, oz, o, w, and a represent Rayleigh scattering, ozone, other gases, water vapor, and aerosols, respectively. 6 Extraterrestrial solar irradiance -- The mean extraterrestrial solar irradiance, H0(8), is taken from the revised Neckel and Labs [1984] data for the wavelength range of 330 to 700 nm. The extraterrestrial solar irradiance corrected for earth-sun distance is given by Gordon et al. [1983] as
2 2 ( JD - 3) F0 ( ) = H 0 ( )1 + ecc cos 365 where ecc is the orbital eccentricity (= ­0.0167) and JD is Julian day of the year. Atmospheric path length -- The slant path length through the atmosphere, M(2), is required for atmospheric transmittance due to attenuation by all constituents. It may be expressed as 1/cos2 for solar zenith angles < 75E, but a correction for the sphericity of the earth-atmosphere system is required at larger zenith angles. Gregg and Carder [1990] used the empirical formulation of Kasten [1966], but we use an updated formulation from Kasten and Young [1989], which is valid at all zenith angles: M ( ) = 1 cos - 0.50572(96.07995 - ) -1.6364 Ozone requires a slightly longer path length for accurate transmittance computations because its dominant concentrations are located in the stratosphere [Paltridge and Platt, 1976]: M oz ( ) = 1.0035 (cos + 0.007)1 / 2
2 Rayleigh scattering -- The Rayleigh total scattering coefficient is taken from Bird and Riordan [1986]: M ( ) Tr ( ) = exp - 4 2 115.6406 - 1.335
where 8 is in :m and M'(2) is the atmospheric path length corrected for atmospheric pressure, M ( ) = M ( ) P P0 P is the atmospheric pressure and P0 is standard atmospheric pressure. The normalized water-leaving 7 radiance (Lwn) algorithm also requires P and will get it from numerical weather models, probably from NMC, according to Dr. Howard Gordon's ATBD-MOD-18. We will take P from the same source. Ozone absorption -- Ozone transmittance is computed via Toz ( ) = exp[- a oz ( ) H oz M oz ( )]
where aoz(8) is the ozone absorption coefficient and Hoz is the ozone scale height. Spectral values of aoz(8) are taken from Inn and Tanaka [1953] and differ slightly from those tabulated by Bird and Riordan [1986] due to the higher spectral resolution here. Hoz should be available as a MODIS product. If not otherwise known, the ozone scale heights can be estimated from the empirical climatological expression of van Heuklon [1979]. Gas and water vapor absorption -- Oxygen is the only atmospheric gas that absorbs significantly in this spectral range. We adopt expressions for transmittance due to oxygen and water vapor absorption from Bird and Riordan [1986]: 1.41 a 0 ( ) M ( ) T0 ( ) = exp- 0.45 [1 + 118.3 a 0 ( ) M ( )] 0.238 a w ( ) WV M ( ) Tw ( ) = exp- 0.45 [1 + 20.07 a w ( ) WV M ( )]
The oxygen and water vapor absorption coefficients (ao and aw, respectively) are derived from transmittance calculations with the 5S Code from Tanre et al. [1990], using the high spectral resolution transmittance observations of Kurucz et al. [1984] to obtain 1-nm resolution. WV is the total precipitable water vapor in cm, which is MODIS product MOD05. Note that the expression for oxygen gas transmittance uses the pressure-corrected path length, M'(2). Aerosol scattering and absorption -- Aerosol concentrations and types vary widely over time and space. Consequently, accurate prediction of their optical thicknesses is difficult. The original Gregg and Carder [1990] model estimated aerosol optical thickness, Ja(8), using the Navy aerosol model [Gathman, 1983], which is parameterized by the local meteorological variables "air-mass type", 24 hr. average wind speed, instantaneous wind speed, and relative humidity. Here, life is simpler because the atmospheric correction procedure for MODIS radiances provides the information necessary to compute Ja(8). 8 First, we write the Angstrom formulation for aerosol optical thickness: a ( ) = -
[Van de Hulst, 1981] where $ is the turbidity coefficient, which is independent of wavelength and (1) represents the aerosol concentration, 8 is wavelength in :m, and " is the Angstrom exponent. We then make a ratio of Eq. 1 at 8 = 412 and 667 nm, take the logarithm, and isolate " on the left to get (412) ln a a (667) = 667 ln 412
Among the atmospheric correction parameters provided in MODIS product MOD37 are Ja(869) and the "epsilon" values, ,(8i,8j), which are defined as: (i , j ) = a (i ) a (i ) p a ( , 0 , i ) a ( j ) a ( j ) p a ( , 0 , j ) where 8i and 8j are any two MODIS wavebands, Ta is the aerosol single-scattering albedo, and pa is the aerosol scattering phase function. For marine or non-absorbing aerosols, the approximation a (412) (412,869) a (667) (667,869)
should be valid [Gordon et al., 1983]. Substitution provides our expression to compute ": ( 2) (412,869 ln (667,869) = 667 ln 412
$ is then calculated via = a (869) 869
" and $ are then used in Eq. 1 to compute Ja(8), and aerosol transmittance is computed by Ta ( ) = exp[- a ( ) M ( )]. 9 The clear-water epsilon product (MOD39, ATBD-MOD-21) flags pixels with highly absorbing aerosols (e.g., Saharan dust). This flag will thus also indicate pixels where the IPAR/ARP products are less accurate, due to the approximation used in Eq. 2. 3.1.2.1.2 Diffuse irradiance ­ Eds(8,0+) Eds(8,0+) is computed via E ds ( ,0 + ) = I r ( ) + I a ( ) + I g ( )
where Ir, Ia, and Ig represent the diffuse components of incident irradiance arising from Rayleigh scattering, aerosol scattering, and multiple ground-air interactions, respectively. Ig is set to zero because multiple sea surface-boundary-layer/atmosphere interactions are rare [Gordon and Castano, 1987]. Rayleigh scattering -- Ir is computed by I r = F0 cos Toz Tu Tw Taa (1 - Tr 0.95 ) 0 .5 (3) (8 dependencies are now suppressed) where Taa represents the transmittance after aerosol absorption (not scattering). All of the other components on the right-hand side of the Eq. 3 are computed in the direct irradiance calculations. Taa is given by Taa = exp[- (1 - a ) a M ( )]
[Justus and Paris, 1985], where Ta is the single-scattering albedo of the aerosol. Ta is computed as a = (-0.0032 AM + 0.972) e 0.000306 RH
where AM is the Navy aerosol model air-mass type and RH is the percent relative humidity. AM ranges from 1 for marine aerosol-dominated conditions to 10 for continental aerosol-dominated conditions. It is assumed to be 1 over the ocean unless the absorbing aerosol flag from MODIS product MOD39 (clearwater epsilon product, ATBD-MOD-21) is set, in which case AM is set to 10. We will get RH from the same source as does the [Lw]N algorithm, which will be the output of numerical weather models, probably from NMC, according to ATBD-MOD-18. Aerosol scattering -- Ia is computed by I a = F0 cos Toz To Tw Taa Tr 1.5 (1 - Tas ) Fa 10 where Tas represents transmittance due to aerosol scattering only and Fa is the forward scattering probability of the aerosol. Tas is computed as Tas = exp[- a a M ( )]
[Justus and Paris, 1985]. Following Bird and Riordan [1986], Fa is computed from the following set of equations: Fa = 1 - 0.5 exp[(B1 + B2 cos ) cos ] B1 = B3 [1.459 + B3 (0.1595 + 0.4129 B3 )] B3 = ln(1- < cos > ) B2 = B3 [0.0783 - B3 (0.3824 + 0.5874 B3 )] is the asymmetry parameter, which is an anisotropy factor for the aerosol scattering phase function as a function of 2 [Tanre et al., 1979]. In this algorithm, is given as a function of the aerosol size distribution and can be parameterized in terms of ": < cos >= -0.1417 + 0.82
For " < 0.0, is set to 0.82, while for " > 1.2, is set to 0.65. This is done so that for low ", typical of maritime conditions, the asymmetry parameter converges to the marine aerosol model of Shettle and Fenn [1979], and for high ", typical of continental conditions, the asymmetry parameter converges to that used by Bird and Riordan [1986]. 3.1.2.2 Calculation of IPAR IPAR is defined as
700 IPAR = 1 - E d ( ,0 ) d hc 400 ( 4) where h is Planck's constant and c is the speed of light. IPAR is calculated from Edd(8,0+) and Eds(8,0+) in two steps. First, the sub-surface irradiances are computed. Then the spectra are added together and integrated over the entire spectrum. The downwelling direct and diffuse irradiances just below the sea surface are given by E dd ( ,0 - ) = E dd ( ,0 + ) (1 - d ) E ds ( ,0 - ) = E ds ( ,0 + ) (1 - s ) 11 where Dd is the direct sea surface reflectance and Ds is the diffuse sea surface reflectance. Total downwelling irradiance just below the sea surface, Ed(8,0­) is simply E d ( ,0 - ) = E dd ( ,0 - ) + E ds ( ,0 - )
3.1.2.2.1 Sea Surface Reflectance Dd and Ds are both composed of two terms, d = dsp + f s = ssp + f
[Koepke, 1984] where Ddsp is the direct specular reflectance, Dssp is the diffuse specular reflectance, and Df is reflectance due to sea foam. In general, the reflectances are functions of 2 and wind speed, but these dependencies have been suppressed for brevity. Df is a function of sea surface roughness, which in turn has been related to wind speed, W [Koepke, 1984]. Using Koepke's [1984] observations, Gregg and Carder [1990] developed the following expressions relating Df to W, which we also use. For W # 4 m s­1, f =0
for 4 < W # 7 m s­1, f = 0.000022 a C D W 2 - 0.00040
C D = 0.00062 + 0.00156W -1
and for W > 7 m s­1, f = (0.000045 a C D - 0.000040)W 2
C D = 0.00049 + 0.000065W
where Da = 1.2 x 103 g m­3 is the density of air and CD is the drag coefficient. The expressions for CD are based on those of Trenberth et al. [1989] and on Koepke's observations that Df = 0 for W # 4 m s­1. Comparing Df calculated by the above equations with Koepke's observations yield a root-mean-square (rms) error of 2.54% for the range 4 to 20 m s­1. By not including foam reflectance, the error in total direct reflectance at 20 m s­1 for a zenith sun was > 52%. By including this formulation, the error was reduced to 1.2%. Foam reflectance is considered isotropic and thus has no dependence on 2. Ddsp is dependent on 2, and for a flat ocean it can be computed directly from Fresnel's law. However, Austin [1974] and Preisendorfer and Mobley [1986] have shown that Ddsp is also dependent on 12 sea state, which can be related to wind speed. Gregg and Carder [1990] developed the following pair of expressions relating Ddsp to 2 and W, which we also use. First, for 2 < 40º or W < 2 m s­1, dsp ( ) = 0.5 sin 2 ( - r ) tan 2 ( - r ) + sin 2 ( + r ) tan 2 ( + r ) where 2 is the solar zenith angle and 2r is the refracted solar zenith angle, which is derived from the expression sin = nw sin r
where nw is the index of refraction for seawater, taken to be 1.341 [Austin, 1974]. Second, for 2 $ 40º and W $ 2 m s­1, dsp = 0.0253 exp[b( - 40)]
b = -0.000714W + 0.0618
which is an empirical formulation derived from Austin's data. This empirical expression is only applied where 2 $ 40º because Fresnel's law is still approximately valid for all wind speeds up to 2 m/s. This formulation produced reflectances within 9.5% rms of the data tabulated by Austin, which, incidentally, also agreed with Preisendorfer and Mobley's ray-tracing calculations to within 10% rms, despite Austin's neglect of multiple reflections. The diffuse specular reflectance Dssp is independent of 2. Assuming a smooth sea and uniform sky, it is given a value of 0.066 [Burt, 1954]. For a wind-roughened surface (W > 4 m/s), Dssp decreases to 0.057 [Burt, 1954]. 3.1.2.2.2 Integration of Ed over wavelength We approximate the integral in Eq. 4 by using a weighted sum at each of the visible MODIS wavelengths. The new formulation for IPAR is IPAR = 1 6 i Ed i ,0 - wEd (i) hc i -1 ( ) where 8i = 412, 443, 488, 531, 551, and 667 nm, and wEd(i) is the weighting function. The Appendix describes the weighting function and its derivation. 13 3.1.2.3 Calculation of ARP The main use of ARP will be as an input to the chlorophyll fluorescence algorithm. Since 90% of the water-leaving radiance is due to scattering in the top attenuation depth [Gordon and McCluney, 1975], we assume that most of the photons fluoresced by chlorophyll which are detected from space also will originate from there. ARP is defined here as
700 z 685 ARP = 400 a ( ) E
0 0 ( , z ) dz d (5) where aN is the phytoplankton absorption coefficient provided by MOD18, E0 is the scalar irradiance, and z685 is calculated as z 685 cos r a w (685) + a (675) (6) aw(685) is the water absorption coefficient at 685 nm, [Carder et al, 1999] taken from Pope and Fry [1997], and aN(675) is taken from the output of the Case 2 chlorophyll algorithm. E0 is E0 ( z) = Ed ( z) µ d ( z) + Eu ( z ) µ u ( z) (7 ) where Ed and Eu are the downwelling and upwelling irradiances and µ d andµ u are the downwelling and upwelling average cosines (the wavelength dependency has been suppressed for brevity). Ed and Eu can be written as E d ( z ) = E d (0 - ) e - K d z , E u ( z ) = E u (0 - ) e - K u z where z is the depth in m and Kd and Ku are the downwelling and upwelling diffuse attenuation coefficients in m­1, both assumed here to be constant over the depth range of interest. For brevity, let's look at ARP at just any given wavelength, eliminating the wavelength integral in Eq. 5. Substituting Eqs. 6 and 7 into Eq. 5 and taking constant terms outside of the depth integral yields ARP = a E d (0 - ) z 685 µd e
0 -Kd z dz + a E u (0 - ) z 685 µu e
0 - Ku z dz 14 Evaluating the two integrals, we get
z 685
0 e - K d z dz = 1 - e - K d z685 Kd z 685 ,
0 e - K u z dz = 1 - e - K u z685 Ku Kd and Ku can be approximated as K d (a + bb ) / µ d , K u (a + bb ) / µ u where a is the total absorption coefficient and bb is the total backscattering coefficient. a is output at the visible MODIS wavelengths by the Case 2 chlorophyll algorithm, and we assume bb