Retrieving Atmospheric Dust Opacity on Mars by Imaging Spectroscopy at Large Angles
RRetrieving Atmospheric Dust Opacity on Mars byImaging Spectroscopy at Large Angles
S. Douté Institut de Planétologie et d’Astrophysique de Grenoble (IPAG), France([email protected] Phone: +33 4 76 51 41 71 Fax +33 4 76 51 41 46)
X. Ceamanos Meteo France CNRM/GMME/VEGEO
T. Apperé Laboratoire AIM CEA-Saclay
Abstract
We propose a new method to retrieve the optical depth of Martian aerosols(AOD) from OMEGA and CRISM hyperspectral imagery at a reference wave-length of 1 µ m. Our method works even if the underlying surface is completelymade of minerals, corresponding to a low contrast between surface and atmo-spheric dust, while being observed at a fixed geometry. Minimizing the effect ofthe surface reflectance properties on the AOD retrieval is the second principal as-set of our method. The method is based on the parametrization of the radiativecoupling between particles and gas determining, with local altimetry, acquisi-tion geometry, and the meteorological situation, the absorption band depth ofgaseous CO . Because the last three factors can be predicted to some extent,we can define a new parameter β that expresses specifically the strength of thegas-aerosols coupling while directly depending on the AOD. Combining estima-tions of β and top of the atmosphere radiance values extracted from the observedspectra within the CO gas band at 2 µ m, we evaluate the AOD and the surfacereflectance by radiative transfer inversion. One should note that practically β can be estimated for a large variety of mineral or icy surfaces with the exceptionof CO ice when its 2 µ m solid band is not sufficiently saturated. Validation Preprint submitted to Elsevier 17th May 2013 a r X i v : . [ a s t r o - ph . E P ] M a y f the proposed method shows that it is reliable if two conditions are fulfilled:(i) the observation conditions provide large incidence or/and emergence angles(ii) the aerosol are vertically well mixed in the atmosphere. Experiments con-ducted on OMEGA nadir looking observations as well as CRISM multi-angularacquisitions with incidence angles higher than 65° in the first case and 33° inthe second case produce very satisfactory results. Finally in a companion paperthe method is applied to monitoring atmospheric dust spring activity at highsouthern latitudes on Mars using OMEGA. Keywords:
Mars; Atmosphere; Radiative Transfer; Aerosols; OMEGA;CRISM
Introduction
Visible and near infrared imaging spectroscopy is a key remote sensing tech-nique to study and monitor the planet Mars. Although its atmosphere is muchfainter than Earth’s, its composition dominated by CO gas implies numerousand strong absorption bands that often overlap with spectral features comingfrom the surface. Furthermore small mineral particles or H O ice clouds of-ten drift over Martian surfaces at various altitudes. These aerosols have alsoa strong, spatially and temporally varying influence on the morphology of theacquired spectra. As a consequence accurate analysis for the study of surfacematerials requires the modeling and the correction of the atmospheric spectraleffects. The first step in this matter consists in retrieving the aerosol opticaldepth (AOD) over the scene.Since 2004 the imaging spectrometer OMEGA aboard Mars Express per-forms nadir-looking and EPF (Emission Phase Function) observations in theVIS (visible) and the SWIR (short wave infrared) (920-5100 nm at 14 to 23 nmspectral resolution) for the study of the surface and the atmosphere of the redplanet. See Clancy and Lee (1991) for the definition of the EPF mode. Thespatial resolution of OMEGA typically varies between 300 and 2000 m/pixeldue to its eccentric orbit. The authors in Vincendon et al. (2007) have de-2eloped a method to quantify the contribution of atmospheric dust in SWIRspectra obtained by OMEGA regardless of the Martian surface composition.Using multi-temporal observations at nadir with significant differences in solarincidence angles, they can infer the AOD and retrieve the surface reflectancespectra free of aerosol contribution. However, this method relies on the veryrestrictive assumption that the atmosphere opacity remains approximately con-stant during the time spanned by the employed acquisitions.In Vincendon et al. (2008) the same authors map the AOD in the SWIRabove the south seasonal cap of Mars from mid-spring to early summer. Thismapping is based on the assumption that the reflectance in the 2.64 µ m satu-rated absorption band of the CO ice at the surface is mainly due to the lightscattered by aerosols above most places of the seasonal cap. In this case, onegeometry is sufficient for the AOD retrieval. Nonetheless, this method is re-stricted to the area of CO deposits that are not significantly contaminated bydust nor water ice.The Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) isa hyperspectral imager on the Mars Reconnaissance Orbiter (MRO) spacecraft.In targeted mode, a gimbaled Optical Sensor Unit (OSU) removes most along-track motion and scans a region of interest that is mapped at full spatial andspectral resolution (18 or 36 m/pixel, 362-3920 nm at 6.55 nm/channel). Inthe targeted mode, ten additional, spatially binned images (180 m/pixel) aretaken over the same region before and after the main image at 10 emergenceangles ranging from -70° to 70°. They provide the so-called Emission PhaseFunction (EPF) for the site of interest that is intended for atmospheric studyand correction of atmospheric effects.Regarding the atmospheric correction of CRISM data, McGuire and 14co-authors (2009) adapted and improved the so-called volcano-scan technique(Langevin et al., 2006). This method removes the CO gas absorption bandsof any spectrum of interest after division by a scaled reference spectrum (i.e.the ratio between the atmospheric transmission at the summit and the base ofOlympus Mons evaluated on a Martian sol when the amounts of ice and dust3erosols were minimal). This simple technique works reasonably well for sur-faces spectrally dominated by minerals, water ice, and sometimes CO2 ice butdoes not correct for aerosol effects .In McGuire and 23 co-authors (2008), a DISORT-based model retrievedthe dust and ice AOD, the surface pressure and temperature from previousexperiment products as well as the acquisition geometry, and the measured I/Fspectrum as inputs. Then, a surface Lambertian albedo spectrum is computedas the output. However this algorithm does not retrieve the AOD directly fromthe images nor does it take advantage of the EPF measurements of the CRISMtargeted mode.Brown and Wolff (2009) proposed a first attempt in that direction by usingthe DISORT algorithm to model the signal at one wavelength (i.e. 0.696 µ m).They iteratively adjust three parameters (surface albedo, dust and ice opacity)in order to achieve a close fit at five points spread across the EPF curve. Never-theless the method is time consuming and the surface albedo is assumed to beLambertian. It has been proved that this assumption bias the AOD and surfaceestimation (Lyapustin, 1999).In this article, we propose an original method that overcomes the previouslimitations to retrieve the optical depth of the Martian dust from OMEGA orCRISM data at the reference wavelength of one micron. The method is basedon a parametrization of the radiative coupling between aerosol particles and gasthat determines, based on the local altimetry and the meteorological situation,the absorption band depth of gaseous CO . We consider the intensity of theabsorption feature at 2 µ m as a proxy of the AOD, provided that the otherinfluencing factors have been taken into account. Our method relies on radia-tive transfer calculations that assume lambertian properties for the surface eventhough, as demonstrated in this paper, the influence of the latter hypothesisis minimized by considering the effect of the aerosols on the gaseous absorp-tion and not the total signal. When processing OMEGA observations, we arecomplementary to the method of Vincendon et al. (2008) since our approachprocesses pixels occupied by a wider variety of materials - pure mineral or water4ce as well as CO and H O deposits contaminated by a large amount of dust -while being observed at a fixed geometry. When processing CRISM observationswe take full advantage of the top of the atmosphere spectral EPF measured bythe instrument for the retrieval of the AOD.This article is organized as follows. In section 1 we describe the post-processing and the formatting of the sequence of EPF calibrated image cubesaccompanying the high resolution CRISM observation. In section 2, we giveinsights about Martian atmospheric properties and radiative transfer (RT) inthe SWIR range. Furthermore we describe models that calculate the spectralradiance coming from Mars and reaching the sensor at the top of the atmo-sphere. In section 3 we expose the basic assumptions, RT parametrization, andproperties on which relies the method for retrieving the AOD from the images.In addition we describe how we implement the method that is then validatedon synthetic data. Finally, in section 4, we present and discuss results obtainedfor representative observations, one from OMEGA and the others from CRISM.Conclusions are eventually drawn in the last Section.
1. Post-processing and formatting of CRISM EPF observations
Due to the composite nature of CRISM EPF acquisitions, such data needsome post-processing prior to AOD retrieval. First of all, CRISM data aretransformed from apparent I/F units (i.e. the ratio of reflected radiance toincident intensity of sunlight) to reflectance units. A Lambertian surface issupposed and data are divided by the cosine of the solar incidence angle (Murchieand 49 co-authors, 2007).CRISM data are corrected for artifacts caused by non-uniformities of the in-strument, residuals of the radiometric correction or external sources. First, hy-perspectral images are corrected for striping and spiking effects. Corrections forboth distortions have been proposed by Parente (2008). Secondly, the CRISMspectrometer is affected by a common artifact to “push-broom" sensors, the so-called “spectral smile" effect. Smile effects refers to the artifacts originated bythe non-uniformity of the instrument spectral response along the across-track5imension, i.e. , the horizontal axis that corresponds to the data columns. Themitigation of smile effects is crucial since the spectral band corresponding to theabsorption maximum of the CO gas at 2 µ m is particularly affected by spectralsmile. We remember that the proposed AOD retrieval method is based on thisspectral feature. As a matter of fact, the estimation of the AOD might be biasedif smile effects are not addressed. Ceamanos et al. developed a twofold methodthat corrects CRISM observations for smile by mimicking an optimal smile-freespectral response (Ceamanos and Doute, 2010). Nevertheless since uncertaintiesassociated to the method itself might propagate to the AOD retrievals we con-servatively restrict the use of each spatially binned image to the central columnsthat own the best spectral response the so-called “sweet spot” columns. Thisdecimation of the data is not a problem because the EPF sequence providesplenty of points for each eleven geometries.CRISM observations are spatially rearranged to evaluate the atmosphereoptical depth of a given terrain position depending on geometry. First, thecentral image that owns the highest spatial resolution is binned by a factor often to match the spatial resolution of the EPF series. After that, all images areprojected onto the same geographical space. At this point, CRISM observationsshowing a high overlap of the eleven acquisitions can be discriminated. In fact,a good overlap assures the existence of a high number of terrain units thathave been sensed from many geometries. In a single hyperspectral data setthat gathers all EPF images and the central scan (hereafter called CSP, CubeSpectroPhotometric), each super pixel conjugated to a given terrain unit gathersup to eleven spectra depending on the completeness of the overlap. More detailscan be found in (Ceamanos et al., 2013).
2. Mars atmospheric radiative transfer
In the present section we give insights about Martian atmospheric propertiesand radiative transfer (RT) in the SWIR range. Furthermore we describe modelsthat calculate the spectral radiance coming from Mars and reaching the sensorat the top of the atmosphere. 6e consider the Martian atmosphere as a vertically stratified medium witha plan parallel or spherical geometry. We divide the atmosphere into ≈ l is at height z l above the surface. Each layer contains a mixture of gas andsolid particles and is homogeneous in temperature T l [K], pressure P l [mbar]and gaseous composition C l [ppm]. A dozen atmospheric absorption bands due to gaseous CO and marginallyto H O are present in the SWIR. We calculate the transmission function of theatmosphere along the vertical, for each pixel of our spectral images, and as afunction of wavelength.
LBLRTM.
For that purpose, we employ a line-by-line radiative transfer model(LBLRTM) (Clough et al., 1995) fed by the vertical compositional and ther-mal profiles predicted by the European Mars Climate Database (EMCD) For-get et al. (2006) for the dates, locations, and altitudes of the observations.The LBLRTM code takes into account the vertical profile of temperature, pres-sure and gaseous composition and interrogates the spectroscopic database HI-TRAN («high-resolution transmission molecular absorption database», ) to compute for each homogeneous layer l the monochromatic absorption coefficient a νl [m ] ( ∆ ν = 10 − cm − ) of thegaseous mixture for the wavenumber ν . The absorption lines typical of eachmolecular specie are characterized by a «Voight» profile that describes at thesame time the broadening by Doppler effect and by collision. The code LBLRTMfrom the Atmospheric and Environmental Research (AER) company also takesinto account various spectroscopic coupling between species (some of themnot being directly active in the optical sense like N or the inert gases). If u l [m -2 ] is the vertically integrated molecular density in terms of number inlayer l , the vertical monochromatic transmission follows the Beer-Lambert law: T ν ( u l ) = exp ( − a νl u l ) ethod k-correlated coefficients. In order to take into account the fine details ofthe gas lines, the final RT calculation, which will integrate the gaseous as wellas particulate absorption and scattering, should be conducted in principle ac-cording to the maximum spectral resolution of the problem ( ∆ ν = 10 − cm − ).Since we must simulate tens, even hundreds, thousand spectra with numeri-cally intensive codes, such a requirement cannot be satisfied practically. Eachspectral channel k of the sensors that we consider (see Introduction) presents awidth ≈
10 cm -1 and thus contains hundreds of lines that cannot be consideredindividually. Consequently, only a statistical approach like the one proposed byGoody et al. (1989) offers a practical solution of the problem. Starting from theline spectra in the spectral interval [ ν k , ν k +1 ] ( ν k , being the lower bound of thecurrent channel k ), we build the discrete form of the frequency distribution forthe absorption: f k ( a i ) = ∆ νν k +1 − ν k W ( a i , a i + ∆ a i )∆ a i (1) a i , i = 1 , . . . , N is a partition of the interval of values taken by the monochro-matic absorption coefficient a νl in [ ν k , ν k +1 ] , W ( a i , a i +∆ a i ) the number of spec-tral points having a value falling in the bin i of lower bound a i and of width ∆ a i .First, we choose an irregular partition that gives the same number of pointsfor each bin i . The objective is to estimate numerically the frequency f k ( a i ) (i.e.a probability) with bins all containing the same number of individuals ( ≈ , anumber large enough to guaranty a correct and homogeneous accuracy. Becausethere are much more lines of weak intensity than lines of strong intensity, thewidth ∆ a i is a strongly increasing function of i . From f k ( a i ) the discrete formof the cumulative frequency distribution g k ( a ) is derived such as: g k ( a n ) = (cid:80) ni =1 f k ( a i )∆ a i . The inverse function a ( g ) = (cid:0) g k (cid:1) − ( a ) , which is generallydefined in the literature as a «k distribution» and is established in the interval g ∈ [0 , , is the highest absorption value reached by the lower fraction g ofthe population of lines ordered according to their intensity. In practice, the8revious function is calculated for any value of g by spline interpolation of theexperimental curve passing through the points (cid:0) g k ( a n ) , a n (cid:1) .Second, we split the line population in ( N + 1 ) sub-sets each representinga fraction ∆ g i = g i +1 − g i , i = 1 , . . . , N − , ∆ g = g and ∆ g N = 1 − g N ofthe total. A Gaussian quadrature g i , i = 1 , . . . , N of order N=2 m , m ∈ N forexample allows an optimal sampling of the weak and strong absorption values.Then, the function a ( g ) allows to build the corresponding partition but in thespace of absorption values a i = a ( g i ) . We calculate the mean value a i of themonochromatic absorption coefficient in each bin [ a i , a i + ∆ a i ] thanks to thevalues a νl returned by LBLRTM AER.Finally, as in Lacis and Oinas (1991), we assume that the mean transmissionvalue for the instrument channel k through the layer l is: T k ( u l ) = (cid:90) ∞ exp ( − au l ) f k ( a ) da = (cid:90) exp ( − a ( g ) u l ) dg ≈ N +1 (cid:88) i =0 exp( − a l,ki u l )∆ g i (2)Furthermore, we introduce a mean optical depth linked to the sub-sets i =0 , . . . , N + 1 , each representing a fraction ∆ g i of the population of lines sortedby increasing intensity : τ gazl,ki = a l,ki u l . Then T l,kgaz = (cid:80) N +1 i =1 w i exp (cid:16) − τ gazl,ki (cid:17) with w i the set of coefficients associated to the Gaussian partition g i . Spectral transmittivity.
If we now consider the whole stack of atmospheric layers,we must make an additional hypothesis in order to calculate an approximationof the total vertical transmission. It is based on the fact that the atmosphericcomposition does not change significantly over the heights z spanned by theregion that contributes predominantly to the transmission. Thus, the form ofthe distribution a ( g ) is very similar from one layer to the other even though theabsolute level -which depends on the pressure level for example- varies. Indeed,if an intense absorption line is centered at wavenumber ν for layer l then theline originating from the same transition has a high probability of falling nearlyat the same position for the layer l (cid:48) (cid:54) = l even though its width and intensityare different. In other words, the probability that a strong absorption level9n the layer l at wave number ν is associated to a weak level in the layer l (cid:48) is very low. Statistically such a situation translates to the following rule con-cerning the probability that the absorption level for the wavenumber ν belongsto the sub-set ∆ g i in the layer l and to the sub-group ∆ g j in the layer l (cid:48) : P (cid:16) (∆ g i ) l | (∆ g j ) l (cid:48) (cid:17) ≈ δ ij (Kronecker symbol). The transmission through bothlayers that should be generally written: T k ( u l , u l (cid:48) ) = N (cid:88) i N (cid:88) j P (cid:16) (∆ g i ) l | (∆ g j ) l (cid:48) (cid:17) exp ( − ( a l,ki u l + a jl (cid:48) ,k u l (cid:48) ))(∆ g i ) l (∆ g j ) l (cid:48) (3)can then be approximated by : T k ( u l , u l (cid:48) ) ≈ N (cid:88) i exp (cid:16) − ( a il,k u l + a il (cid:48) ,k u l (cid:48) ) (cid:17) ∆ g i (4)since we use the same quadratures for both layers. This formula can be gen-eralized to the whole stack of layers and the vertical transmission through theatmospheric gases in a multilayered medium is given by: T kgaz = N +1 (cid:88) i ∆ g i exp( − (cid:88) l ¯ τ l,kgaz i ) (5)Because the complete calculation is very time consuming, it is only performedfor a limited ( (cid:46) ) number of reference points representative of “regions” slicedaccording to regular bins in latitude lat , longitude long and altitude h . Threetransmission spectra are generated for the maximum, mean and minimum alti-tude of a given region. Then the transmission spectrum of all the pixels belong-ing to the region is interpolated from the triplet depending on their individualaltitude to give T kgaz ( h, lat, long ) Optical properties of Martian mineral aerosols are still largely undocumentedeven though recent studies have improved our understanding (Korablev et al.,2005). We favor the single scattering albedo, optical depth spectral shape andphase function retrieved in the near-infrared by Vincendon et al. (2008) and10olff et al. (2009). In the first case,a Henyey-Greensten phase function withone parameter is used. This is a coarse parametrization but it is relevant forthe phase angle range spanned by the data set of nadir OMEGA observationsthat we consider. In the second case, because it has been proved that theHenyey-Greenstein function is somewhat inaccurate to recreate the effect ofaerosols for the CRISM angular ranges, the radiative properties are describedwith more details. A single scattering albedo, optical depth spectral shape, andphase function based on cylindrical particles are considered. The compositionof aerosols has been shown to be remarkably homogeneous, e.g. between Roverlanding sites. Changes of radiative properties (Vincendon et al., 2009) mainlyresult from differences in the mean of the aerosols size distribution (Wolff et al.,2003) because the latter results from a balance between sedimentation and liftingmechanisms, which in turn depend on turbulence, density, etc. Additionally thedegree of hydration for the mineral particles might change with latitude andseason . Fortunately, for the typical hydration values reached on Mars, thespectral effects are only noticeable between 2.7 and 3.5 µ m (Pommerol et al.,2009), a range that we thus ignore. Water ice crystals may also be frequentlypresent but are not considered in this paper because they may influence theCO gas 2 µ m absorption feature in a specific manner. Because the atmosphere isusually well mixed, we take a vertical distribution of optical depth at a referencewavelength ( k : 1 µ m) that is exponentially decreasing with height such that: τ l,k aer = τ k aer (exp( − z l /H scale ) − exp ( − z l +1 /H scale )) / (1 − exp( − zmax/H scale )) τ k aero = (cid:80) Ll =1 τ l,k aer is the column integrated opacity of the aerosols at the refer-ence channel k , H scale is the scale height of the distribution, z l the height oflayer l lower interface and zmax = 100 km. The scale height H scale and theoptical depth τ k aero remain the only free parameters concerning the atmosphere. In general, photons undergo absorption and multiple scattering events wheninteracting both with gas and aerosols. One should note that we neglect Rayleigh11cattering due to atmospheric molecules because, according to the formula ofBodhaine et al. (1999), that will lead for a mean ground pressure of 6 mbarto an optical depth of τ ray ≈ . − - far lower than the contribution we canexpect for the aerosols in any case. The value of the monochromatic aerosoloptical depth (at ∆ ν = 10 − cm − ) does not vary significantly within a channel k whereas the counterpart for the gas can vary by several orders of magnitude.As a consequence calculating at channel k the transfer of radiation in the atmo-sphere requires to couple gas and aerosols properties at each atmospheric level l and for each sub-set i = 0 , . . . , N + 1 of the gas line population. As indicatedin Table 1, the coupling procedure leads to N + 1 sets of optical depth, singlescattering albedo and phase function that characterize the layers. These prop-erties, as well as the acquisition geometry ( θ i , θ e , φ e ), are introduced one by onein a RT engine (e.g. DISORT, Stamnes et al. (1988)) in order to calculate foreach level z the partial radiance field I ki ( z, θ i , θ e , φ e ) corresponding to each sub-set lines+aerosols. The surface is characterized by its bidirectional reflectancedistribution function (BRDF). The total radiance field is the mean value of thepartial fields weighted by the fractions ∆ g i : I k ( z, θ i , θ e , φ e ) = N +1 (cid:88) i ∆ g i I ki ( z, θ i , θ e , φ e ) (6)
3. Method for retrieving the optical depth
The proposed method is based on a parametrization of the radiative couplingbetween aerosol particles and gas that determines, with local altimetry andthe meteorological situation, the absorption band depth of gaseous CO . Thecoupling depends on (i) the acquisition geometry (ii) the type, abundance andvertical distribution of particles (iii) the bidirectional reflectance factor of thesurface (BRF). For each spectrum of an image, we compare the depth of the2 µ m absorption band of gaseous CO that we estimate on the one hand fromthe observed spectrum and on the other hand from a calculated transmissionspectrum through the atmospheric gases alone. This leads to a relevant new12arameter that directly depends on τ k aer . Combining the latter with the radiancefactor R k obs within the 2 µ m band, we evaluate by RT inversion the AOD τ k aer and the reflectance factor of the surface. Our method is based on two main assumptions:1. First we assume that the surface reflects the solar and atmospheric ra-diation isotropically and, consequently, that the ground-atmosphere interface ischaracterized by a single scalar quantity (the normal albedo A ksurf ) that dependson wavelength.2. Second we parametrize the radiative coupling between aerosols and gasby assuming that the latter contribute to the signal as a simple multiplicativefilter. In this way, the top of the atmosphere (TOA) radiance is written suchthat I k ( θ i , θ e , φ e ) ≈ T kgaz ( h, lat, long ) (cid:15) k ( θ i ,θ e ,φ e ,τ k aer ,H scale ,A ksurf ) I ksurf + aer ( θ i , θ e , φ e ) (7)the effect of the acquisition geometry, the aerosols, and the surface Lamber-tian albedo A ksurf being taken into account by scaling the aerosol free verticaltransmission T kgaz ( h, lat, long ) -calculated according to 2.1- by an exponent (cid:15) . I ksurf + aer ( θ i , θ e , φ e ) is the calculated TOA radiance but without considering theatmospheric gases in the radiative transfer. Factor (cid:15) k can be further decomposedinto two terms: (cid:15) k ( θ i , θ e , φ e , τ k aer , H scale , A ksurf ) = ψ ( ν ) β k ( θ i , θ e , φ e , τ k aer , H scale , A ksurf ) (8)with the airmass ν = θ i ) + θ e ) . Factor ψ is purely geometric and allowsa quick and simplified calculation of the free gaseous transmission along theacquisition pathlength (incident-emergent) using the vertical free transmission: 13 (cid:48) kgaz = (cid:88) i a i exp( − (cid:88) l ¯ τ l,kgaz i .ν ) ≈ ( T kgaz ) ψ ( ν ) (9)the approximate formula being ψ ( ν ) = p ( p + p ν + p ν + p ν ) . In orderto retrieve the values of the factors p to p , we fit the experimental cloud ofpoints obtained by plotting log (cid:16) T (cid:48) kgaz (cid:17) / log ( T kgaz ) as a function of the airmass ν for a set of representative atmospheric vertical profiles. For that purpose thequantities T (cid:48) kgaz and T kgaz are respectively the slant and vertical versions of thefree transmission through the gas given by the LBLRTM code. Note that thedispersion around the mean curve is very limited. . Factor p accounts for thedifference of spectral resolution and radiometric calibration between OMEGAand CRISM and its retrieval is explained at the end of section 3.5. All p-factorvalues are given in Table 2.Factor β k on the other hand expresses the aerosol effect on gaseous ab-sorption and depends on the considered channel k . In practice this spectraldependability can be easily managed. Indeed, numerical experiments and theanalysis of real data (see following sections) show that this exponent dependsprincipally on the gaseous absorption intensity T kgaz . If (i) β k can be estimatedexperimentally ( ˆ β k ) for one or for several geometries for each pixel (super-pixel)of an OMEGA (CRISM) observation and (ii) if one assumes a value for A k surf ( k is the spectral band number at 2 µ m) then ( τ k aer ) can be derived. Indeed,function β k ( θ i , θ e , φ e , τ k aer , H scale , A k surf ) is invertible provided that the scaleheight H scale of the aerosols is known. In this matter, several studies such asVincendon et al. (2008) suggest that H scale = 11 km is a good guess, but insection 4 we prove that H scale = 8 km may be more appropriate as it leads tobetter results. Our optical depth retrieval system uses a multi-dimensional look-up table(LUT) of β values arranged according to discrete combinations of acquisitiongeometries ( θ i , θ e , φ e ) and of physical parameters ( τ k aer , H scale , A k surf ) , see Table3. For each combination, the value of β is calculated on the basis of three TOA14pectra. T kref is a reference free vertical transmission corresponding to a giventime and location on the surface of Mars. It has been chosen in order to cover arange of absorption values (all spectral bands considered) typical of the southernhemisphere in spring. The spectra I k and I ksurf + aer in units of radiance aregenerated successively by the code DISORT with and without considering thegaseous absorption in the atmosphere. A series of k-correlated coefficients ¯ τ l,kref i intervenes in the calculation of T kref and I k . Finally the LUT is built accordingto the formula: β ref ( θ i , θ e , φ e , τ k aer , H scale , A k surf ) = k (cid:48)(cid:48) (cid:88) k = k (cid:48) f k β k = k (cid:48)(cid:48) (cid:88) k = k (cid:48) f k ln( I k I ksurf + aer ) / ln( T (cid:48) kref ) (10)with f k = 1 /K . The channels k = k (cid:48) , . . . , k (cid:48)(cid:48) ( K = k (cid:48)(cid:48) − k (cid:48) + 1 ) encompassthe 2 µ m CO gas absorption band such as T kref < . . For DISORT calcula-tions, we consider the atmosphere as a plan parallel vertically stratified mediumwhich is a sufficient approximation to calculate the TOA radiance wheneverincidence and emergence directions are ≈ o above the horizon. Several exper-imental scatter-plots T kref - β k built using the reference and also other sets ofk-correlated coefficients show that there exists a linear relationship between thetwo quantities with an excellent correlation factor ( (cid:38) . see figure 1) : β k = a ( ν, τ k aer , H scale , A k surf ) + b ( ν, τ k aer , H scale , A k surf ) T kref (11)We notice that the coefficients ( a, b ) of the regression depend weakly on theparameters τ k aer , H scale , A k surf . Because the intensity of the CO gas 2 µ m absorption feature is diverselyaccessible from the spectra depending on the nature of the surface materialspresent in the pixel, we need to consider different strategies for estimating thefactor β . the albedo of the surfaces that we consider rarely exceeds 0.6 ineral surfaces: method 1. The factor β can be readily estimated for everyOMEGA or CRISM spectrum that does not show any features, such as the H Oor CO ice signature, superimposed on the 2 µ m CO gas absorption band. Inthis way, this absorption band is completely due to the atmosphere. A similarformula to the one used in the volcano scan technique (McGuire and 14 co-authors, 2009) is used replacing the Olympus reference transmission spectrumby ( T kgaz ) ψ ( ν ) such that ˆ β k = αψ ( ν ) ln (cid:32) TkgazTk gaz (cid:33) α = ln (cid:18) R kobs R k obs (cid:19) + 0 . (cid:18) R k obs R k obs (cid:19) (12)with k the index of a channel falling at the shortwave extremity of the 2 µ m band wing. The scatter-plot T kgaz - ˆ β k , k = k (cid:48) , . . . , k (cid:48)(cid:48) built for each pixelof a real OMEGA observation confirms the existence of a linear relationshipbetween the two quantities with a generally moderate dispersion (see Figure1). The regression coefficients (ˆ a, ˆ b ) are then calculated and will allow us toestimate exponent β in the same absorption range than the reference : ˆ β ref = ˆ a + ˆ b k (cid:48)(cid:48) (cid:88) k = k (cid:48) f k T kref (13)Note that we use the series of channels k = k (cid:48) , . . . , k (cid:48)(cid:48) for the estimationof ˆ β ref allowing us to be less sensitive to noise and channel aging than usingthe single channel k . In order to perform the inversion process, we take as aninitial value for the lambertian albedo of the surface A k surf = R k obs ( T k gaz ) ψ ( ν )ˆ βk . Water ice surfaces: method 2.
A specific procedure is necessary for any spec-trum marked by the signature of water ice with the possible presence of dustbut without spectral contamination by CO ice. A numerical optimization isperformed regarding a cost function that depends on β and that expresses thequality of the gaseous absorption correction on the spectra. The quality criteriawe choose is the band shape of solid H O at 2 µ m that must show a unique local16inimum and a simple curvature. Mathematically this criteria is equivalent tothe second derivative of the spectrum that must be close to zero for channels [ k (cid:48) , k (cid:48)(cid:48) ] covering the H O band and slightly beyond. After all the cost functionis written: ξ ( β ) = k (cid:48)(cid:48) (cid:88) k = k (cid:48) (cid:32) d dk R kobs ( T kgaz ) ψ ( ν ) β (cid:33) (14)Searching for the global minimum of this function leads to the estimation ofthe exponent β : ˆ β = arg min β ξ ( β ) . The normalization of ˆ β up to the referencelevel of absorption uses the regression line: ˆ β ref = ˆ b ( k (cid:48)(cid:48) (cid:88) k = k (cid:48) f k T kref − T k gaz ) + ˆ β Once again, in order to realize the inversion we take as initial value for thelambertian albedo A k surf = R k obs ( T k gaz ) ψ ( ν )ˆ β . Carbon dioxide ice (CO ). Currently there is no simple strategy to estimate in apixelwise manner the β coefficient from spectra presenting the signature of solidCO because the latter cannot be distinguished from its gaseous counterpart at2 µ m. For this kind of surface, our method to estimate τ k aer is ineffective.If, on the one hand, CO ice proves to be sufficiently pure with grain sizesmore than 100 µ m in diameter (this conditions is always respected in practicefor the seasonal deposits) then the method proposed by (Vincendon et al., 2008)is the only solution. It is based on the assumption that the reflectance level at2.63 µ m (channel k ) is nearly equal to zero because of the absorption saturationof the ice. Consequently, in this case, the reflectance factor R k obs represents thepath radiance of the atmosphere (mostly due to the aerosols at channel k )and is a bijective, easily invertible function of τ k aer . In section 3.5 we will seehow exponent β can be calculated from τ k aer , computed using (Vincendon et al.,2008), and R k obs for correcting the gaseous absorption on the spectra (section3.5). 17f, on the other hand, CO ice proves to be substantially contaminated bydust or water ice, then the only conceivable way to estimate β could be astatistical procedure as it is proposed in Luo et al. (2010).Sulfates present a broad band similar to water ice around 2 µ m and areoccasionally spectrally dominant in OMEGA or CRISM observations. A futurerequirement is finding a way to deconvolve the CO gaseous features from thesulfates 1.94 strong absorption in order to evaluate β .Hereafter, we simplify the writing of the parameters by removing subscript ref from the notations. The quantities β and ˆ β are as if they were wavelengthindependent because they are normalized up to the reference level of absorption. We now examine the sensitivity of function β regarding itsdifferent parameters θ i , θ e , φ e , τ k aer , H scale , A k surf . This study is summarized asfollows. The exponent β decreases with the airmass ν especially as the aerosolopacity is on the rise. This is the direct consequence of an intensifying couplingeffect between gas and aerosols, the latter reducing the effective pathlengthof the photons in the atmosphere by scattering and thus the probability tobe absorbed by CO gas. The exponent evolution reflects such a decrease ofrelative absorption band intensity. The coupling is higher when the sensor looksin the solar direction (90 (cid:46) φ e (cid:46) °) than when its looks in the opposite halfspace (0 (cid:46) φ e (cid:46) °) . The higher the scale height characterizing the verticaldistribution of the aerosols, the more aerosols influence atmospheric absorptionespecially since their opacity is high. The surface reflectance factor also controlsthe absorption by CO gas because of occurrence of multiple scattering betweenthe two media. A bright surface promotes a higher number of roundtrips forthe photons which then have a higher probability to be absorbed: exponent β increases with Lambertian albedo A k surf . Finally, because the sensitivity of β with τ k aer ( ∂β∂τ ) decreases with the airmass, we can expect significant estimationerrors for τ k aer below a certain threshold. Indeed the ground pressure providedby the Martian climate database (Forget et al., 2006) is only indicative of the18eal pressure at the moment of the observation. We must stress that the groundpressure determines T kgaz then ˆ β and, for this reason, an uncertainty ∆ p surf willpropagate as an uncertainty ∆ β on the coupling factor, then as an uncertaintyon the optical depth ∆ τ k aer = ∂τ∂β ∆ β . We quantify this effect in the validationsection 3.6.The experimental behavior of ˆ β , evaluated above mineral surfaces basedon real data with respect to airmass ν and τ k aer , is also full of information.We distinguish two situations according to whether we work with OMEGA orCRISM observations. OMEGA data.
In the OMEGA case, we consider large populations of pixelsextracted from several hyperspectral images. The principal variability of theaerosol opacity across the scene comes from altitude changes -sometimes severalkilometers in elevation- that modify the local atmospheric height and thus thecolumn integrated density of the aerosols. The horizontal heterogeneity of theparticle density at constant altitude generally implies a more moderate vari-ability. Hence, by selecting pixels belonging to slices spanning no more than acouple hundred meters in elevation but observed under various geometries, wecan draw a scatter-plots ( ν, ˆ β ) such as that presented in Fig. 2. Set apart a no-ticeable internal dispersion linked to factors such as those previously mentionedor spatial variations of Lambertian surface albedo, the points correspondingto the OMEGA observation 1880_1 are organized according to monotonouslydecreasing curves of ˆ β over a wide range of increasing ν values. The curvesfaithfully follow the theoretical scatter plot that can be built thanks to the nu-merical reference LUT for H scale ≈
11 km and A k surf ≈ . . The best fit givesan estimation of the mean optical depth τ k aer ≈ CRISM data.
The sensitivity study is conducted on CRISM EPF sequenceFRT82EB. This observation corresponds to an ice-free cratered area of the southpole of Mars ( lat =-83 ◦ ). The image geometry corresponds to the conditions ofthe high latitudes with θ i =74.5 ◦ . Due to the sun-synchronous orbit of CRISM,there are two modes in azimuthal angle, φ e (cid:39) ◦ and φ e (cid:39) ◦ . Regarding the19opography, FRT82EB presents mild slopes and a constant altitude. The initialalbedo of the surface can be approximated by the average value of the spectralband at 1 micron. In this case, the A k surf of FRT82EB is equal to 0.3.Data are processed by the data pipeline described in Section 1, except thefiltering of the columns affected by the “spectral smile” effect. Thanks to thegood overlap of the high-resolution image and of the EPF, we can define up to850 super pixels that are represented by more than 6 geometries.The factor β can be readily estimated for every CRISM spectrum that doesnot show CO ice signature superimposed on the 2 µ m CO gas absorptionband. Once a series of spectra corresponding to a super pixel is treated, weobtain β as a function of up to eleven geometries. Figure 3 summarizes theresults for the image FRT82EB by plotting β as the function of the airmassdistinctively for the two azimuths explored by the observation (black and greycrosses).One can see how the azimuthal angle has an impact on the gas-aerosol cou-pling due to aerosol phase function strong anisotropy (Wolff et al., 2009). Inaddition, the model curves of β that provide the best matching with the exper-imental results are shown by plain lines. In fact, the best match correspondsto a value of τ k aer around 0.5, respectively 8 km for H scale . Curves for othertargeted observations are shown in Fig. 10. We now describe the practical implementation of the method that allows thecomputation of surface albedo and aerosol opacity maps from a given observa-tion. The implementation comes in two flavors depending on the origin of thedata.
OMEGA.
The processing of a nadir OMEGA image is conducted on a pixelwisebasis after segmenting the image according to an automatic detection of thesurface components Schmidt et al. (2007). For the pixels that do not presentthe signatures of solid CO the following iterative algorithm is carried out:20 nput : a spectrum R kobs of rank p in the flattened image, the correspondingvertical transmission T kgaz calculated as explained in section 2.1, the acquisitiongeometry ( θ i , θ e , φ e ) and a global estimation of the scale height H scale
1. estimating the exponent β (a) using method 1 (section 3.3) for a mineral surface(b) using method 2 for a water ice surface without CO ice contamination2. initialization i = 0 , Lambertian albedo (cid:16) A k surf (cid:17) = R k obs ( T k gaz ) ψ ( ν )ˆ β
3. DO4. building a 1D curve τ k aer as a function of β by quadrilinear interpolationof the LUT described in section 3.3 for a given set ( θ i , θ e , φ e ) , A k surf = (cid:16) A k surf (cid:17) i and H scale .5. linear interpolation of the previous curve at ˆ β to calculate (cid:0) τ k aer (cid:1) i +1
6. building a 1D curve A k surf as a function of R k T OA by quadrilinear inter-polation of the LUT for a given set of ( θ i , θ e , φ e ) and τ k aer = (cid:0) τ k aer (cid:1) i +1
7. linear interpolation of the previous curve at R k obs ( T k gaz ) ψ ( ν )ˆ β in order to cal-culate (cid:16) A k surf (cid:17) i +1 i = i + 1
9. UNTIL (cid:18)(cid:12)(cid:12)(cid:12)(cid:12)(cid:16) A k surf (cid:17) i − (cid:16) A k surf (cid:17) i − (cid:12)(cid:12)(cid:12)(cid:12) ≤ .
01 AND (cid:12)(cid:12)(cid:12)(cid:0) τ k aer (cid:1) i − (cid:0) τ k aer (cid:1) i − (cid:12)(cid:12)(cid:12) ≤ . (cid:19) OR i ≥ imax The algorithm usually converges after a dozen iterations and leads to a completesolution once A k surf is calculated from R k obs and τ k aer based on a dedicated LUT.If the initial estimation ˆ β or, to a lesser extent, the estimation of (cid:16) A k surf (cid:17) ispoor then the convergence is not reached and there is no solution.For the pixels presenting a signature of solid CO pure enough (accordingto the criteria established by Vincendon et al. (2008)) we have seen in section3.3 that it is possible to obtain τ k aer by their independent method. The problemis then laid down differently: we must calculate the exponent β from τ k aer and R k obs that will be used for correcting the gaseous absorption on the spectra. Thealgorithm is also iterative in this case :21 nput : a spectrum R kobs of rank p in the flattened image, the correspondingvertical transmission T kgaz extracted from the ATM cube at the same rank,the acquisition geometry ( θ i , θ e , φ e ) and a global estimation of the scale height H scale
1. estimation of τ k aer using the method by Vincendon et al. (2008)2. initialization i = 0 , Lambertian albedo (cid:16) A k surf (cid:17) = 0 .
3. DO4. building a 1D curve β as a function of τ k aer by quadrilinear interpolationof the LUT described in section 3.3 for a given set ( θ i , θ e , φ e ) , A k surf = (cid:16) A k surf (cid:17) i and H scale .5. linear interpolation of the previous curve at τ k aer to calculate (cid:16) ˆ β (cid:17) i +1
6. calculation of C obs = R k obs ( T k gaz ) ψ ( ν ) ( ˆ β ) i +1
7. building a 1D curve A k surf as a function R k T OA by quadrilinear interpo-lation of the LUT for a given set ( θ i , θ e , φ ) and τ k aer
8. linear interpolation of the previous curve at C obs in order to calculate (cid:16) A k surf (cid:17) i +1 i = i + 1
10. UNTIL (cid:18)(cid:12)(cid:12)(cid:12)(cid:12)(cid:16) A k surf (cid:17) i − (cid:16) A k surf (cid:17) i − (cid:12)(cid:12)(cid:12)(cid:12) ≤ .
01 AND (cid:12)(cid:12)(cid:12)(cid:12)(cid:16) ˆ β (cid:17) i − (cid:16) ˆ β (cid:17) i − (cid:12)(cid:12)(cid:12)(cid:12) ≤ . (cid:19) OR i ≥ imax CRISM.
The processing is applied to the single hyperspectral data set thatgathers the central scan and the EPF images, i.e. the CSP cube. We choosean area of the image that offers the best compromise between its size and itsspectral homogeneity. In addition, this area should not present the signatures ofsolid CO . All the super-pixels that fall in the area are blended into one vectorof parameters that needs to be explained as a whole by the model provided thesolution for parameters τ k aer and R k obs is found. The following iterative algorithmis carried out. 22 nput : the collection of valid spectra R kobs , the corresponding vertical transmis-sions T kgaz extracted from the ATM cube at the same coordinates, the acquisitiongeometries ( θ i , θ e , φ e ) and a global estimation of the scale height H scale
1. for each valid pixel p , estimating the exponent β (a) using method 1 (section 3.3) for a mineral surface(b) using method 2 for a water ice surface without CO ice contamination2. initialization i = 0 , mean Lambertian albedo for the selected area (cid:16) A k surf (cid:17)
3. DO4. building a matrix of β values, each element β pq corresponding to thecoupling factor that is predicted by the model for the spectrum of index p assuming that the atmospheric opacity is (cid:0) τ k aer (cid:1) q , q = 1 , . . . , Q , a samplingof the range of possible variation for this parameter. The value β pq iscalculated by quadrilinear interpolation of the LUT described in section3.3 for a given set ( θ i , θ e , φ e ) , A k surf = (cid:16) A k surf (cid:17) i , (cid:0) τ k aer (cid:1) q , and H scale .5. finding the global minimum of the objective function χ ( τ k aer ) = (cid:80) p =1 ,...,P (cid:16) ˆ β p − β pq (cid:17) to estimate (cid:0) τ k aer (cid:1) i +1
6. for each valid pixel p , building a 1D curve A k surf as a function of R k T OA by quadrilinear interpolation of the LUT for a given set of ( θ i , θ e , φ e ) and τ k aer = (cid:0) τ k aer (cid:1) i +1 . Linear interpolation of the previous curve at R k obs ( T k gaz ) ψ ( ν )ˆ β in order to calculate (cid:16) A k surf (cid:17) i +1
7. calculation of the mean Lambertian albedo (cid:16) A k surf (cid:17) i +1 for the selectedarea8. i = i + 1
9. UNTIL (cid:18)(cid:12)(cid:12)(cid:12)(cid:12)(cid:16) A k surf (cid:17) i − (cid:16) A k surf (cid:17) i − (cid:12)(cid:12)(cid:12)(cid:12) ≤ .
01 AND (cid:12)(cid:12)(cid:12)(cid:0) τ k aer (cid:1) i − (cid:0) τ k aer (cid:1) i − (cid:12)(cid:12)(cid:12) ≤ . (cid:19) OR i ≥ imax Note that the value of factor p appearing in Table 2 is adjusted for the OMEGA(respectively CRISM) dataset by maximizing the overall convergence rate achievedby algorithms 3.5a (respectively 3.5c) for a selection of observations. Note thatthis optimization problem is satisfactorily convex and thus its solution satis-23actorily constrained, since the convergence of the iterative inversion is verysensitive to the value of p . The validation is realized by inverting synthetic data, that is to say, real-istic TOA spectra I k calculated with DISORT (section 2.3) to which we addzero-mean Gaussian noise simulated using a meaningful covariance matrix, i.e.calculated from an estimation of OMEGA noise. In the reference simulation, weset the properties of the surface A k surf = 0 . (minerals) and of the atmosphere τ k aer = 0 . , H scale = 11 km, and initial pressure p surf . In a second simula-tion, for the calculation of the spectrum I k , we apply a pressure deviation of ∆ p surf = ± Pa typical of the Martian meteorological variability (baroclineactivity, (Forget et al., 1999)) but the factor β is estimated with a transmissionspectrum T kgaz calculated according to the initial pressure. In a third simu-lation, we perturb the initial distribution profile of the aerosols (exponential)by increasing the opacity of one given layer l, l = 1 , . . . in succession by anadditional amount ∆ τ k aer = 0 . . Thus, the total AOD becomes 0.7 with theperturbed profile and there are as many runs as the layer number. In each run,the geometry varies according to the values given in Table 3. In agreement withthe sensitivity study (section 3.4), Figure 4 shows that our method estimates τ k aer with an accuracy comparable to the usual methods in the reference caseprovided that the airmass exceeds ν ≈
3. Since we do not have access to thereal ground pressure but to a value predicted by a MGCM, we expose ourselvesto potentially important errors if ν (cid:46) but lower than 10% beyond ( (cid:46) ν (cid:46) ).Figure 5 illustrates the fact that a perturbation of opacity does not have thesame impact on factor β and thus on the estimation τ k aer according to its heightabove the ground. A low lying aerosol layer is not detectable by our methodbecause it does not have much influence on the mean pathlength of the photonsinto the gas. By contrast, a detached layer high in the atmosphere superimposedon the standard profile has a strong influence on the coupling. Thus, such alayer could skew our estimation because the latter does not take into account24uch a local perturbation. In conclusion, we state that our method is reliable iftwo conditions are full filled: (i) the observation conditions provide an airmassthat is large enough, i.e. ν (cid:38) . (ii) the aerosol particles are well verticallymixed caused by a vigorous convection.
4. Results
Figure 6 shows in a synthetic way all the products obtained by the chainof operations described above for the observation OMEGA
ORB1880_1 that isrepresentative of many. The chosen image covers a large geographical area ofMars at high and medium southern latitudes in spring. The bright extendedpart appearing in the maps R k obs and A k surf represents the seasonal deposits offrozen CO as well as the permanent cap that is, at that time, buried by thelatter. The black coded area indicate the surfaces (among then a significantfraction of the “cryptic” region - where CO very much contaminated by dustexists - for which no method for estimating τ k aer is currently available or indicatethe pixels for which the iterative inversion has not worked (for example near thelimb).First cross validation of the two methods used in order to produce the map τ k aer - ours and and the one by Vincendon et al. (2008)- can be performed asfollow. We plot along-track profiles of optical depth through the maps and weexamine to which extent the segments corresponding to the mineral surfaces (inblack) and those corresponding to the frozen surfaces (in grey) are well con-nected in figure 7. The right-side red segment with high dispersion correspondsto the “cryptic” region. The quality of the connection seems quite good for ORB1880_1 and for other observations that we treated.Second the problem of evaluating the dust scale height, which is an input ofour retrieval algorithms, is investigated more thoroughly. Our region of interestis the southern high latitudes in spring for which dust activity is monitoredin the companion paper. In this case there are two possibilities for constrain-ing the dust scale height provided that it can be considered spatially constant25or one given observation. First its value (along with factor p ) can be ad-justed by maximizing the overall convergence rate achieved by algorithm 3.5awhen analyzing each image of a selection of OMEGA observations. In all cases H scale =11 km turns out to be the optimal value. Note that this optimizationproblem is satisfactorily convex, and thus its solution satisfactorily constrained,since the convergence of the iterative inversion is very sensitive to the value of( p , H scale ). It is important to keep in mind that the previous scale heightexpresses the rate of optical depth change with height. Second large scale geo-graphical variations of τ k aer can be interpreted as due, at first order, to changesof the atmospheric column linked with changes of pixel altitude. Indeed if weassume that the intrinsic optical properties (such as the extinction cross-section σ ) and the scale height of dust particles are both geographically and verticallyconstant and that their density number N normalized to altitude 0 is alsogeographically constant, then τ k aer ( h ) = N σH exp ( − h/H ) = τ exp ( − h/H ) .We are restricted to OMEGA observations that cover an area with a range ofaltitudes large enough so that the experimental cloud of points obtained byplotting τ k aer as a function of the pixel altitude h can be satisfactorily fitted bythe model. Figure 8 shows the result for four global observations. The data canbe satisfactorily explained with values of H varying from 6 to 11 kilometers,most often in the lower end of this interval, in agreement with Vincendon andLangevin (2010). We note that the “spatially” derived scale height H does notcorrespond to the “radiatively” derived scale height H scale ( H (cid:46) H scale ) withthe exception of observation ORB1849_1. If we trust the method a possibleway to explain this discrepancy is to suppose that N must show a overall sys-tematic trend of decrease with altitude. In a sense this parameter expressesthe intensity of dust loading in the atmosphere regardless of the length of theatmospheric column. We hypothesize that the systematic trend is linked withthe decrease capability of the atmosphere to lift dust with the decrease of at-mospheric density while the dispersion around the main trend is linked with themeteorological activity. If we assume also an exponential altitude profile witha scale height of H (cid:48) scale for N , the effective “spatially” derived scale height for26 k aer is H ≈ H scale H (cid:48) scale / H scale + H (cid:48) scale ≈ H scale = H (cid:48) scale =11 km, con-sistent with what we found. Referring to Figure 8, it is interesting to note thatthe observation with the highest spatially derived scale height corresponds to apeak of dust activity around L S ≈ ° also observed by Pancam during MY27(Lemmon and Athena Science Team, 2006). This global dust event erases theoverall systematic trend of N decrease with altitude.Because it is impossible to constrain H scale for each OMEGA observation wechoose to take a fixed value H scale =11 km. We have determined that possibleexcursions of H scale by ± image and an RGB com-position made by extracting the TOA martian reflectivity measured by OMEGAat three wavelengths (0.7070, 0.5508, and 0.4760 µ m) from the correspondingvisible image (Figure.9). The composition is stretched so as to make visiblethe dust clouds as yellowish hues against mineral surfaces while the seasonaldeposits appear completed saturated. Excellent qualitative agreement can benoted between the two products down to the details. In the companion paperother successful examples are provided. In the CRISM case, six observations of the landing site of the Mars Explo-ration Rover “Spirit” (i.e. West of Columbia Hills), acquired with varied ge-ometry and atmospheric conditions (see Table 4), have been chosen to test ourmethod (hereafter identified as the “ β method”). For that purpose, these scenesare of special interest since the AOD values were also retrieved from the samedata by Wolff et. al. (personal communication) and from PanCam synchronousmeasurements at the ground looking towards zenith (Johnson et al., 2006). Theprevious procedures are identified as the “W method” and the “P method” re-spectively. Special attention is paid by focusing on areas of the observations withaltitude variations not exceeding one hundred meters and with negligible slope27t the scale of the pixel size and above. Following the procedure described insection 1, the eleven CRISM images of each FRT observation are combined andbinned at about 300 meters per pixel. Then, the algorithm adapted to CRISMfor the simultaneous retrieval of surface albedo and aerosol opacity is applied onthe CSP cube (section 1). The operation leads to the results tabulated in Table4 along with those obtained by the “W” and “P” methods. In addition Figure 10illustrates to which extent our best model matches the experimental ˆ β versusairmass curves for each observation. Similarly to the OMEGA case, airmass ν ≈ represents a threshold above which the fit becomes very satisfactory withthe exception of the observation FRT95B8. One should also note that the largerthe difference between the two azimuthal modes is, the more constrained thebest solution likely becomes because the two branches of the ˆ β versus airmasscurves are increasingly separated.Once the content of aerosols is known, CRISM TOA radiances can be cor-rected in order to retrieve surface BRF. The correction of remotely sensed im-ages for atmospheric effects is however not straightforward due to the anisotropicscattering properties of the atmospheric aerosols and the materials at the sur-face. Traditional inversion algorithms are based on reductionist hypothesis thatassumes that the surface is lambertian (McGuire and 23 co-authors, 2008). Al-though this assumption largely simplifies the inverse problem, it critically cor-rupts the angular shape of the retrieved BRF since solid surfaces are hardlyisotropic (Lyapustin, 1999). Recently, we have proposed an original inversionmethod to overcome these limitations when treating CRISM multi-angle obser-vations Ceamanos et al. (2013). This inversion algorithm is based on a TOAradiance model that depends on the Green’s function of the atmosphere and asemi-analytical expression of the surface BRF. In this way, robust and fast in-versions of the model on CRISM TOA radiance curves are performed accountingfor the anisotropy of the aerosols and the surface. The root mean square errorachieved by the model when reproducing the CRISM TOA reflectance factorcurves is indicated in Table 4 for each observation and for each AOD input. ANaN value indicates that the inversion was not successful, i.e. no valid surface28RF could be retrieved with the proposed aerosol properties and atmosphericopacity.By examining Table 4 we conclude that our method gives results generallyin good agreement with the values retrieved by Wolff et al. with the notableexception of observations 812F and 95B8. In the first case, the “W” methodoverestimates τ k aer since, with such an opacity, the inversion algorithm has diffi-culties to produce a physically meaningful surface reflectance. At contrary ourvalue is in agreement with the one estimated from the PanCam measurements.In the second case, the ˆ β versus airmass curves are not well fitted by our modelimplying that the estimated value for τ k aer by the “ β method” may not be ac-curate. One should also note that the root mean square error achieved by theTOA reflectance factor modeling using our value is systematically lower thanwhen using the value given by the “W method” except for observation 334D.We put forward the hypothesis that the AOD estimate by the “ β method” isless biased even though the latter also uses the lambertian surface hypothesisin the radiative transfer calculations. As regards to the previous hypothesis,the relatively weak dependence of the main parameter ˆ β from which the AODis retrieved is a possible explanation. Indeed the radiative coupling betweenaerosols and gas is mostly expressed in the atmospheric additive term and, toa lesser extent, the terms including multiple scattering of photons between sur-face and atmosphere. In the first case at least the anisotropy of the surface hasno influence. By contrast, methods that directly invert a lambertian surface-atmosphere RT model on the TOA radiance in the continuum of the spectra,such as the “W method” lead to an AOD estimate that is more sensitive tothe surface bidirectionnal reflectance. We explain the moderate discrepancies ofresults that can be observed between the “P method” and the “ β method” if weconsider that the first one is more sensitive to the low lying layers of aerosolswhereas the second one is more sensitive to aerosols layers situated at ≈ onclusions In this article we propose a method to retrieve the optical depth τ k aer of Mar-tian aerosols from OMEGA and CRISM hyperspectral imagery. The methodis based on parametrization of the radiative coupling between particles and gasdetermining, with local altimetry, acquisition geometry, and the meteorologicalsituation, the absorption band depth of gaseous CO . The coupling depends on(i) the acquisition geometry (ii) the type, abundance and vertical distribution ofparticles, and (iii) the surface albedo A ksurf . For each spectrum of an image, wecompare the depth of the 2 µ m absorption band of gaseous CO between (i) theobserved spectrum and (ii) the corresponding transmission spectrum throughthe atmospheric gases alone. The latter is calculated with a line-by-line RTmodel fed by the vertical compositional and thermal profiles predicted by theEuropean Mars Climate Database for the dates, locations, and altitudes of theobservations. Thus we define a significant new parameter β that expresses thestrength of the gas-aerosols coupling while directly depending on τ k aer . Com-bining β and the radiance value R k obs within the 2 µ m band, we evaluate τ k aer and A ksurf by radiative transfer inversion and provided that the radiative prop-erties of the aerosols are known from previous studies and that an independentestimation of the scale height of the aerosols is available. One should note thatpractically β can be estimated for a large variety of mineral or icy surfaceswith the exception of CO ice when its 2 microns solid band is not sufficientlysaturated.The validation of the method was performed both with synthetic and realdata. It shows that our method is reliable if two conditions are fulfilled: (i)the observation conditions provide large incidence or/and emergence angles (ii)the aerosol are vertically well mixed in the atmosphere. A low lying aerosollayer is not detectable by our method because it does not have much influenceon the mean pathlength of the photons into the gas. By contrast, a detachedlayer high in the atmosphere superimposed on the standard profile has a stronginfluence on the coupling. Our method works even if the underlying surface is30ompletely made of minerals, corresponding to a low contrast between surfaceand atmospheric dust, while being observed at a fixed geometry. This is the firstprincipal asset of our method. Minimizing the effect of the surface Lambertianassumption on the AOD retrieval is the second principal asset of our method.Experiments conducted on OMEGA nadir looking observations -as well asCRISM EPF- produce very satisfactory results. With OMEGA, we note a goodcoherency between our approach and the one of Vincendon et al. (2008). Thedomain of airmass ( ν (cid:38) . ) for which our method is reliable implies that it isintended for analyzing high latitude OMEGA observations with sufficiently highsolar zenith angles ( (cid:38) τ k aer retrievals for solar incidenceangle down to 33° extending the applicability of the method to non polar re-gions. Finally we should note that our method was applied for the first timeextensively on a series of OMEGA observations in order to map the atmosphericdust opacity at high southern latitudes from early to late spring of Martian Year27. This study is presented in the companion paper. Acknowledgments
This work was done within the framework of the Vahiné project fundedby the “Agence Nationale de la Recherche” (ANR) and the “Centre d’EtudesSpatiales” (CNES).
References
Bodhaine, B. A., Wood, N. B., Dutton, E. G., Slusser, J. R., 1999. On rayleigh optical depthcalculations. Journal of Atmospheric and Oceanic Technology 16 (11), 1854–1861.URL http://journals.ametsoc.org/doi/abs/10.1175/1520-0426%281999%29016%3C1854%3AORODC%3E2.0.CO%3B2
Brown, A. J., Wolff, M. J., Mar. 2009. Atmospheric Modeling of the Martian Polar Regions: OneMars Year of CRISM EPF Observations of the South Pole. In: Lunar and Planetary InstituteScience Conference Abstracts. Vol. 40 of Lunar and Planetary Institute Science ConferenceAbstracts. pp. 1675–+.Ceamanos, X., Doute, S., 2010. Spectral smile correction of crism/mro hyperspectral images.Geoscience and Remote Sensing, IEEE Transactions on 48 (11), 3951 –3959. eamanos, X., Douté, S., Fernando, J., Schmidt, F., Pinet, P., Lyapustin, A., 2013. Surfacereflectance of mars observed by crism/mro: 1. multi-angle approach for retrieval of surfacereflectance from crism observations (mars-reco). Jounal of Geophysical Research Planets 118,1–20.Clancy, R. T., Lee, S. W., Sep. 1991. A new look at dust and clouds in the Mars atmosphere -Analysis of emission-phase-function sequences from global Viking IRTM observations. Icarus93, 135–158.Clough, A., S., Iacono, J., M., 1995. Line-by-line calculation of atmospheric fluxes and coolingrates 2. Application to carbon dioxide, ozone, methane, nitrous oxide and the halocarbons.J. Geophys. Res. 100 (9), 16519–16536.Forget, F., Hourdin, F., Fournier, R., Hourdin, andTalagrand, C., O., Collins, M., Lewis, R., S.,Read, andHuot, P. L., J., oct 1999. Improved general circulation models of the Martianatmosphere from the surface to above 80 km. J. Geophys. Res. 104 (.13), 24155–24176.Forget, F., Millour, E., Lebonnois, S., Montabone, L., Dassas, K., Lewis, S. R., Read, P. L.,López-Valverde, M. A., González-Galindo, F., Montmessin, F., Lefèvre, F., Desjean, M.-C.,Huot, J.-P., Feb. 2006. The new Mars climate database. In: Mars Atmosphere Modelling andObservations. pp. 128–+.Goody, R., West, R., Chen, L., Crisp, D., 1989. The correlated-k method for radiation calculationsin nonhomogeneous atmospheres. Journal of Quantitative Spectroscopy and Radiative Transfer42 (6), 539 – 550.URL Johnson, J. R., Grundy, W. M., Lemmon, M. T., III, J. F. B., Johnson, M. J., Deen, R., Arvidson,R. E., Farrand, W. H., Guinness, E., Hayes, A. G., Herkenhoff, K. E., IV, F. S., Soderblom, J.,Squyres, S., 2006. Spectrophotometric properties of materials observed by Pancam on the MarsExploration Rovers: 2. Opportunity. Journal of Geophysical Research 111 (E12S16,).Korablev, O., Moroz, V., Petrova, E., Rodin, A., 2005. Optical properties of dust and the opacityof the martian atmosphere. Advances in Space Research 35 (1), 21–30.URL
Lacis, A. A., Oinas, V., May 1991. A description of the correlated-k distribution method formodelling nongray gaseous absorption, thermal emission, and multiple scattering in verticallyinhomogeneous atmospheres. J. Geophys. Res. 96, 9027–9064.Langevin, Y., Doute, S., Vincendon, M., Poulet, F., Bibring, J.-P., Gondet, B., Schmitt, B.,Forget, F., Aug. 2006. No signature of clear co2 ice from the /‘cryptic/’ regions in mars’ southseasonal polar cap. Nature 442 (7104), 790–792.URL http://dx.doi.org/10.1038/nature05012
Lemmon, M. T., Athena Science Team, Mar. 2006. Mars Exploration Rover Atmospheric Imaging:Dust Storms, Dust Devils, Dust Everywhere. In: Mackwell, S., Stansbery, E. (Eds.), 37thAnnual Lunar and Planetary Science Conference. Vol. 37 of Lunar and Planetary Inst.Technical Report. p. 2181.Luo, B., Ceamanos, X., Doute, S., Chanussot, J., Aug. 2010. Martian aerosol abundanceestimation based on unmixing of hyperspectral imagery. In: Hyperspectral Image and SignalProcessing: Evolution in Remote Sensing, 2010. WHISPERS ’10. IEEE, pp. 1 –4.Lyapustin, A. I., Feb. 1999. Atmospheric and geometrical effectson land surface albedo. Journal ofGeophysical Research 104 (D4), 4127–4143.McGuire, P. C., 14 co-authors, Jun. 2009. An improvement to the volcano-scan algorithm foratmospheric correction of CRISM and OMEGA spectral data. Planetary and Space Science 57,809–815.McGuire, P. C., 23 co-authors, Dec. 2008. MRO/CRISM Retrieval of Surface Lambert Albedos forMultispectral Mapping of Mars With DISORT-Based Radiative Transfer Modeling: Phase 1Using Historical Climatology for Temperatures, Aerosol Optical Depths, and AtmosphericPressures. IEEE Transactions on Geoscience and Remote Sensing 46, 4020–4040. urchie, S., 49 co-authors, May 2007. Compact Reconnaissance Imaging Spectrometer for Mars(CRISM) on Mars Reconnaissance Orbiter (MRO). Journal of Geophysical Research (Planets)112, 5–+.Parente, M., Mar. 2008. A New Approach to Denoising CRISM Images. In: Lunar and PlanetaryInstitute Science Conference Abstracts. Vol. 39 of Lunar and Planetary Inst. Technical Report.pp. 2528–+.Pommerol, A., Schmitt, B., Beck, P., Brissaud, O., Nov. 2009. Water sorption on martian regolithanalogs: Thermodynamics and near-infrared reflectance spectroscopy. Icarus 204, 114–136.Schmidt, F., Doute, S., Schmitt, B., 2007. Wavanglet: An efficient supervised classifier forhyperspectral images. Geoscience and Remote Sensing, IEEE Transactions 45 (5), 1374–1385.Stamnes, K., Tsay, S.-C., Jayaweera, K., Wiscombe, W., 1988. Numerically stable algorithm fordiscrete-ordinate-method radiative transfer in multiple scattering and emitting layered media.Applied Optics 27, 2502–2509.Vincendon, M., Langevin, Y., Jun. 2010. A spherical Monte-Carlo model of aerosols: Validationand first applications to Mars and Titan. Icarus 207, 923–931.Vincendon, M., Langevin, Y., Poulet, F., Bibring, J.-P., Gondet, B., Jul. 2007. Recovery of surfacereflectance spectra and evaluation of the optical depth of aerosols in the near-IR using a MonteCarlo approach: Application to the OMEGA observations of high-latitude regions of Mars.Journal of Geophysical Research (Planets) 112 (E11), 8–+.Vincendon, M., Langevin, Y., Poulet, F., Bibring, J.-P., Gondet, B., Jouglet, D., OMEGA Team,Aug. 2008. Dust aerosols above the south polar cap of Mars as seen by OMEGA. Icarus 196,488–505.Vincendon, M., Langevin, Y., Poulet, F., Pommerol, A., Wolff, M., Bibring, J., Gondet, B.,Jouglet, D., Apr. 2009. Yearly and seasonal variations of low albedo surfaces on Mars in theOMEGA/MEx dataset: Constraints on aerosols properties and dust deposits. Icarus 200,395–405.Wolff, J., M., Clancy, T., R., sep 2003. Constraints on the size of Martian aerosols from ThermalEmission Spectrometer observations. Journal of Geophysical Research (Planets) (E9), 1–1.Wolff, M. J., Smith, M. D., Clancy, R. T., Arvidson, R., Kahre, M., Seelos, F., Murchie, S.,Savijärvi, H., Jun. 2009. Wavelength dependence of dust aerosol single scattering albedo asobserved by the Compact Reconnaissance Imaging Spectrometer. Journal of GeophysicalResearch (Planets) 114, 0–+. ω gas = 0 ω aer Optical depth τ gas i , ∆ g i , i = 0 , . . . , N + 1 τ aer Phase function Υ gas = 0 Υ aer τ i = τ gas i + τ aer , ω i = ω aer / (1 + τ gas i τ aer ) , i = 0 , . . . , N + 1 Υ = Υ aer Table 1: Sum up of the physical quantities characterizing the “elementary” processes of the ra-diative transfer taking place in each homogeneous atmospheric layer. We omit in the notationthe indices l and k for readability. p p p p OMEGA 1.2 0.66312127 0.44429186 -0.024559039 0.00068174262CRISM 0.95 0.66312127 0.44429186 -0.024559039 0.00068174262
Table 2: Value of the coefficients allowing to calculate ψ ( ν ) . i (°) 0 20 40 50 60 70 72 75 78 80 83 θ e (°) 0 2 20 40 50 55 60 70 φ e (°) 0 30 60 90 120 150 180 τ k aer A k surf H scale (km) 4 6 8 11 14 17 20 Table 3: Series of values taken by the incident, emergent and azimuth angles as well as by thephysical parameters used for the construction of the reference LUT τ k aer “W method” (1 µ m) 0.35 1.74 1.92 0.55 0.35 0.33RMSE w/model 0.877E-02 2.39E-02 1.58E-02 1.12E-02 2.47E-02 2.21E-02 τ k aer “P method” (1 µ m) 0.871 1.32 0.93 0.66 0.48 0.31RMSE w/model NaN 2.02E-02 1.05E-02 1.42E-02 1.55E-02 2.38E-02 τ k aer “ β method” (1 µ m) 0.46 1.41 1.00 0.36 0.46 0.38RMSE w/model 1.10E-02 2.13E-02 1.10E-02 0.877E-02 1.65E-02 1.56E-02 Table 4: Optical depth at 1 µ m retrieved by three different methods for our selection of obser-vations. The Root Mean Square Error (RMSE) refers to the adequacy of a non Lambertiansurface-atmosphere radiative transfer model fed by the latter optical depth in reproducing theTOA reflectance factors measured by CRISM. igure 1: Experimental scatter-plots between T kref and β k for different geometries and at-mospheric conditions. The points represented by the * symbols and fitted by regression linesof the same grey level, come from synthetic data. The points represented by the + symbolsand fitted by a single regression line come from the estimation of factor β k for three differentobserved spectra, each of them corresponding to a specific grey level (see section 3.3) .50.60.70.80.91 3 4 5 6 7 8 9 E s ti m a t e d f ac t o r b e t a Airmass
ORB1880 Model
Figure 2: Scatter plot ( airmass ν, ˆ β ) for the observation 1880_1 related to the modeled ( ν, β )for H scale ≈ km, τ k aer ≈ . and A k surf ≈ . . igure 3: Experimental β curves from the CRISM image FRT82EB according to airmass . Themodel curves that provide the best match are also plotted as solid lines. Note the dispersionof β values likely linked with the “spectral smile” affecting the spatially binned images, i.e.the clusters of points in the scatter plot. A e r o s o l op ac it y Airmass dP=-15PadP=0PadP=+15Pa
Figure 4: Influence of the meteorological variability on the estimation of τ k aer depending onairmass ν . A e r o s o l op ac it y Airmass dtau=0.1, z=0km z=1.81 z=5.19 z=8.26 z=20.57 z=31.44 z=88
Figure 5: Influence on the estimation τ k aer of perturbing the reference vertical distributionprofile of the aerosols by a local perturbation of opacity 0.1 as a function of the airmass ν . igure 6: Parameter maps for the OMEGA observation . From top to bottom R k obs , ˆ β , τ k aer and A k surf . E s ti m a t e d ae r o s o l op ac it y Line number
Figure 7: A MEX along track-profile of τ k aer extracted from the aerosol optical depth map1880_1. The method presented in this paper (factor β ) produces the black segment whileusing the method by Vincendon et al. (2008) leads to the grey segment igure 8: Scatter plot of τ k aer as a function of altitude h for the OMEGA observations 1765_1,1786_1, 1849_1, and 1880_1 ( dots). The continuous lines are the result of applying aregression to the data with an exponential τ exp ( − h/H ) , the scale height H varying aroundthe value that insures the best fit (in black). igure 9: Global OMEGA observation . The composition on the left displays an RGBcomposition of the TOA martian reflectivity in the visible which is stretched so as to revealdust in the atmosphere as yellowish hues. The composition on the right displays the AerosolOptical Depth map at 1 micron. igure 10: Experimental β curves from CRISM selected EPF sequences as a function ofairmass plotted distinctively for the two azimuths explored by the observation (black andgrey crosses). In each case, the model curves that provide the best match are also plotted assolid lines.curves from CRISM selected EPF sequences as a function ofairmass plotted distinctively for the two azimuths explored by the observation (black andgrey crosses). In each case, the model curves that provide the best match are also plotted assolid lines.