Interactive comment on “ The impact of surface reflectance variability on total column differential absorption LiDAR measurements of atmospheric CO 2 ” by J . P

face region have been instigated. The factors derived are the relative difference between the surface reflectance variability from the model perspective using MODIS data, and the variability from a spaceborne DIAL perspective using 30 m resolution Landsat data with overlapping footprint viewing geometry. During the orbit simulations these scaling factors are used on the MODIS variability calculations to obtain the correct surface reflectance variability across the footprint pairs. This approach provides the correct magnitude of the surface reflectance for a nadir pointing DIAL system owing to the use of the MODIS BRDF data, as well as the correct reflectance variability observed by the overlapping DIAL footprints through the application of the derived scaling factors.


Introduction
Carbon dioxide (CO 2 ) is a naturally and anthropogenically occurring greenhouse gas in the Earth's atmosphere whose concentration has increased over the last 200 years by approximately 30% (Keeling et al., 1995).Measurements ob-Correspondence to: P. S. Monks (p.s.monks@le.ac.uk) tained from ice cores provide an historical account of atmospheric CO 2 concentrations as far back as 650,000 years, and via comparison with modern in-situ measurements it has been shown that the atmospheric CO 2 concentration is now over 90 ppmv higher than it has ever been in that time (Siegenthaler et al., 2005).The recent and sudden growth observed has been attributed to anthropogenic activities such as the burning of fossil fuels, cement production and land-use change (IPCC, 2007).
Predicting the effect of increased CO 2 levels on the climate requires improved understanding of the terrestrial processes that control carbon fluxes, as well as any potential feedback mechanisms on these processes which may develop as a consequence of climate change (Sarmiento and Gruber, 2002).Tighter constraints on transport models which at present are only capable of calculating fluxes on continental scales are required (Gibert et al., 2004).To better constrain the present models more observations are needed on a global scale to supplement the sparsely located in-situ measurement networks.Satellite remote sensing has been shown to provide the denser sampling required by inversion modeling and an improvement in flux estimation may be achieved if the precision of total column measurements averaged over monthly time scales on an 8 • x10 • footprint is less than 2.5 ppmv (Rayner and O'Brien, 2001).There are a number of satellite instruments at present measuring CO 2 from space including the SCIAMACHY instrument on board the ESA ENVISAT satellite (Bovensmann et al., 1999), the AIRS instrument onboard the NASA AQUA satellite (Gautier et al., 2003), the IASI instrument onboard the METOP-A satellite (Chalon et al., 2001) and the GOSAT instrument launched by JAXA in January 2009 (Hamazaki et al., 2005).Whilst some of these systems may have shown the potential to meet the required spatial coverage and precision to improve flux estimates on continental scales (Barkley et al., 2006;Buchwitz et al., 2006), they also have limitations.In particular passive measurements which use the short wave infrared are limited by their reliance on solar illumination which restricts their latitudinal coverage, and thermal infrared systems are not sufficiently sensitive near the ground where the largest fluxes occur.Furthermore, passive remote sensing systems involve retrieval complexities which suffer from aerosol contamination and radiation path length uncertainties.
Active remote sensing techniques such as total column differential absorption LiDAR (TC-DIAL) offer an alternative without some of these sources of uncertainty and has been the focus of recent papers and feasibility studies (Amediek et al., 2009;Dufour and Breon, 2003;Ehret et al., 2008;Ehret and Kiemle, 2005;Gibert et al., 2004;Loth et al., 2005;Ismail et al., 2004).Many TC-DIAL advantages over current passive remote sensing systems are a product of its use of pulsed (active) radiation.The path length of the returned pulses may be determined via the time between their transmission and reception, measurements can be made equally well at all latitudes both day and night, and by using a differential retrieval method, aerosol and thin cloud interference is avoided.These advantages make TC-DIAL a flexible CO 2 measuring technique capable of diurnal sampling with potential for near-surface sensitivity and therefore an attractive supplement to the current measurement systems (Loth et al., 2005).At present there are no CO 2 TC-DIAL instruments in space, however breadboard demonstrators and aircraft campaigns are being carried out to develop the technology for space application, e.g.(Amediek et al., 2009).
The principle of TC-DIAL is the differential retrieval of scattered intensities from two laser transmissions of similar wavelengths known as the on and off-line pulses.The on-line pulse frequency is very strongly absorbed by CO 2 whereas the off-line pulse frequency is very weakly absorbed by CO 2 (Figure 1).Although the two transmissions have very different absorption cross sections they are closely situated in frequency space to minimise any differences in their aerosol and water vapour interactions.The two transmissions which make up a single measurement are closely located but do not completely overlap Two spectral lines have been identified as particularly appropriate for CO 2 DIAL using temperature sensitivity, strength and interference from other atmospheric species as selection criteria (Loth et al., 2005).The two lines identified were 4875.748957 and 6367.223459cm −1 .The signal returned to the detector with which a CO 2 measurement is made is reflected from aerosols, cloud and the Earth's surface.For TC-DIAL only the latter applies by definition, however instruments which operate by this principle may also be able to use backscatter from clouds to retrieve above cloud CO 2 if the clouds are optically thick.The measurement pair consisting of the on-line and off-line returned intensities is used in a simple retrieval equation to calculate the volume mixing ratio (VMR) (Ehret et al., 2008) (equations 1 to 3).
where N CO2 is the CO 2 VMR, ∆τ is the differential total column optical depth, S on and S of f are the received on and off line intensities respectively, k is the Boltzmanns constant, T is the atmospheric temperature, P is the atmospheric pressure, ρ w is the atmospheric water vapour concentration, σ on and σ of f are the absorption cross sections for the online and offline pulses respectively and l is the orbital altitude.
The received intensities S on and S of f provide the only information available to the TC-DIAL retrieval regarding the probed volume of atmosphere.Any modulation in the ratio of their intensities which does not arise from atmospheric CO 2 absorption creates errors in the differential optical depth.One such source of modulation is signal sensitivity to variations in the Earth's surface reflectance caused by imperfect surface footprint co-location (see Figure 1).
The co-location ambiguity present in each sounding is a consequence of two independent factors.Firstly, laser pointing jitter within the spacecraft generates random offsets in the relative positions of the on and off-line footprints leading to an overlap geometry which varies with each sounding.Secondly, a time delay between the on and off-line transmissions results in a predictable forward footprint separation owing to the satellites orbital velocity.The delay is designed to avoid significant S on and S of f modulation which would arise through the simultaneous detection of multiple pulses.
This paper investigates the uncertainties incurred by the surface reflectance variability on TC-DIAL retrievals over agriculture in the US, England, France and the Czech Republic.A LiDAR model configured with the instrument and satellite specifications is given in Table 2.
The methods and results are compared to those from work by Amediek et al. (2009), in which the problem of surface reflectance variability on TC-DIAL retrievals has also been investigated (Amediek et al., 2009).A fundamental difference between these papers is the method by which the TC-DIAL viewing geometry is achieved.Amediek et al. (2009) tackled a viewing geometry under sampling problem by upscaling narrow swath, high resolution airborne LiDAR measurements to full scale TC-DIAL surface footprints.In contrast, this work uses lower resolution spaceborne data which is downscaled to simulate the TC-DIAL viewing geometry.Both methods use assumptions and some level of interpolation to obtain realistic results and therefore their agreement is desirable.

Methodology
The Leicester LiDAR model (LLM) is a computer simulation of physical processes associated with pulsed TC-DIAL systems, including interactions with aerosols, absorbing media and the Earth's surface.The LLM is linked with a low Earth orbit simulator for nadir pointing laser systems which defines the position in space and time of each TC-DIAL sounding.

Model Surface
The surface interaction is modeled using a bi-directional reflectance distribution function (BRDF) which is modified to account for the hot spot effect and DIAL viewing geometry (Equation 4).The parameterisation for the function comes from the 500 m resolution MODIS BRDF data product MCD43A1.5 (Strahler et al., 1999), which for the present study required a co-addition of three data sets for cloud clearing purposes.The dates used were the 5 th , 15 th and 25 th of May 2007.The MODIS MCD43 product has been shown to have an RMS error of 0.0130 (equivalent to <5%) from continuous field observations of surface albedo at a number of measurement locations (Salomon et al., 2006).
The BRDF model used is a sum of three parameters (Roujean et al., 1992), f iso , f vol and f geo , which combine to give surface reflectivity R. Two of the parameters have associated kernels, the volumetric kernel K vol and the geometric kernel K geo which provide the BRDF's directional components (Wanner et al., 1995).
The volumetric kernel (equation 5) is modified to account for the hot-spot effect (Maignan et al., 2003), a characteristic of DIAL transmission and viewing geometry (Equation 7). where, θ s is the illumination azimuth angle, θ v is the viewing azimuth angle, and φ is the illumination zenith angle.
The geometric kernel is given as where, The result of applying the BRDF model is a global surface reflectance map at 500 m resolution in units of sr −1 .An example is given over Europe in Figure 2.

Resolution scaling
The resolution of the MCD43 product imposes limitations on the models representation of the surface reflectance variability on the scale of TC-DIAL viewing geometry.Variations in reflectance are expected to be observed over less than 50 m in this study, yet the MODIS data product provides an average figure over 500 m.To compensate a scaling factor is derived from higher resolution data to mediate the gap between the MODIS resolution and that of the reflectance variability observed by a spaceborne TC-DIAL instrument.
Landsat 7 radiance data at 30 m spatial resolution are used to provide the required higher resolution surface reflectance scenes.The Landsat 7 instrument is a sun-synchronous scanning radiometer with an equatorial overpass of between 10 and 11 am local time.A consequence of the Landsat orbit is the solar zenith angle of the recorded scenes is sufficiently low to produce shadowing in the data from tall objects on the surface.To avoid biases in deriving the magnitude of surface reflectivity from this data the scenes selected must avoid all significant shadow producing objects such as tall buildings, forests and cliffs.Rural agriculture on flat land is a surface type which is sufficiently shadow free for the purpose of this study as obstructions such as trees and hedge rows are sufficiently short and sparse to be of substantial significance.
Careful vetting of the Landsat radiance data removes areas of settlements, forestry and any cloud obscuration within the scenes.Aerosol contamination is corrected for using the computer software ENVI.
The TC-DIAL viewing geometry of the on and off-line transmissions creates two closely overlapping circles of laser light on the surface.The energy across the footprints is distributed as symmetrical 2 dimensional Gaussian distributions.This non-linearity in the spatial energy density of the laser light allows variations in surface reflectance at any point within the footprints to cause CO 2 retrieval ambiguities.The scaling factors for the various agricultural scenes are derived with this viewing geometry in mind.
Each 120 km 2 Landsat agricultural scene is divided up into 57,600 individual 500 m MODIS pixel sized areas (MPAs).The 256 complete 30 m Landsat samples that occur within each of the MPAs are averaged to simulate the reflectance observed by individual 500 m resolution measurements.The variability in the surface reflectance at 500 m resolution is simulated by calculating the relative difference of two vertically adjacent MPAs.
An arrangement of Landsat samples within each MPA pair is selected by a Landsat arrangement template (LAT).The template is set out to best match the footprint areas of the TC-DIAL viewing geometry.Two LAT's of equal dimensions are defined to simulate the on and off-line transmissions with one positioned a single Landsat pixel higher than the other to simulate a 30 m footprint separation distance.Each LAT is multiplied by a 2 dimensional Gaussian weighting function to account for the power distribution across the footprints.The Gaussian weighted LAT pairs are moved around within each MPA pair to all 324 possible positions.The reflectance difference within each pair is recorded to produce the higher resolution data (Figure 3).This process is carried out for all 28,800 MPA pairs and the difference between the variability observed by the MPAs and that of the higher resolution LAT's within the MPAs provides the scaling factors.
The scaling factors derived using the Landsat resolution are limited to having a minimum footprint separation distance of 30 m, however it is expected that the separation distance of a TC-DIAL system will be no more than approximately 10 m.Two methods of compensation have been considered to reduce the footprint overlap from 30 m to 10 m within the confines of the available resolution.
The first method uses data from (Amediek et al., 2009), whose TC-DIAL aircraft campaign included a study into the effect of footprint separation distance on the observed surface reflectance.The study involved upscaling approaches to increase the footprint size from their airborne TROPOLEX system to that of a spaceborne system.They indicated a factor of approximately 3.5 for the difference observed in surface reflectance between footprints 10 m apart and 30 m apart.Their spaceborne system had an assumed footprint diameter of 100 m and consequently the variability reduction factor is only directly applicable to this diameter.Unfortunately the 30 m scaling factors generated by the LLM use a footprint diameter of 150 m to coincide better with the Landsat resolution.In the absence of further information it is assumed that upscaling the diameter of the footprints from 100 to 150 m will have relatively little effect on the variability reduction factor as both footprints will scale equally with the same separation distance.The reduction factor of 3.5 is therefore used to generate a series of scaling factors which are presented alongside the standard 30 m scaling factors in table 1.
The second method uses linear interpolation between Landsat pixels to allow one of the LATs to shift by only a 1/3 of a pixel as opposed to a whole pixel.This allows the second footprint to be displaced by only 10 m from the first.In a similar sense to the first method this too is carried out in the absence of higher resolution surface reflectance data, and the generated scaling factors are presented alongside the standard 30 m scaling factors in Table 1.

Reflectance variability
To accommodate the resolution scaling factor F , the surface reflectance for each on and off-line footprint is calculated using equations 14, 15 and 16.
S on and S of f are unitless surface reflectance values used by the LLM to define the quantity of returning laser light from the surface to the detector.R on and R of f are the simulated laser intensities reflected back from the surface in units of sr −1 .A and l are the receiver mirror area and orbital altitude respectively which form the solid angle between the surface and the satellite.F is a scaling factor defined to account for the TC-DIAL viewing geometry and its spatial resolution difference to the surface reflectance data, and V is the surface reflectance variability between vertically adjacent MODIS data points in units of sr −1 .

Model Spectroscopy
To determine the consequence of surface reflectance variability on a CO 2 VMR retrieval, correct spectroscopy and realistic atmospheric profiles are required.For this purpose the LLM uses spectral line shape calculations based on HITRAN 2008 line centres and associated parameters to define the position and behaviour of the relevant spectral lines (Lafferty et al., 2009).
A convolution algorithm constrained by these parameters is applied to the lines with a pressure and temperature dependant Voigt line shape (Mitchell, 1971).The result is vertically varying spectral lines incorporating collision and pressure broadening and pressure shift that are used to define vertically resolved absorption cross sections (Equation 17). where Ref is the spectral lines half width half maximum for a reference atmospheric temperature and pressure.For HI-TRAN and therefore the LLM, these points are 101,325 Pa and 273 K. υ is the wavenumber of the laser transmission, υ 0 is the centre lines wave number (4875.74896cm −1 ), P is the atmospheric pressure and T is the atmospheric temperature.
The atmospheric profiles used for the Voigt line shape calculations are latitude dependant MIPAS reference atmospheres (Remedios et al., 2007).
The vertically resolved absorption cross sections define the atmospheric absorption as a function of altitude for the laser transmissions.The LLM simulates each atmospheric level as a discrete entity, calculating the transmission through each level via the Beer-Lambert law (equation 22).
where I i is the intensity of light having passed through the atmospheric level with index i, I i−1 is the transmitted light from the previous level, n is the atmospheric concentration of CO 2 in molecules cm −3 , σ is the absorption cross section in cm 2 , and l is the length of each atmospheric level in cm.

Model configuration
The LLM is configured using parameters partly based on the system definition of the proposed A-SCOPE mission (A-SCOPE, 2008).These parameters are presented in Table 2.
The vast majority of the MODIS pixels within the regions given in Table 1 are sampled to avoid the possibility of a scaling factor bias.This is achieved by the orbit simulator being run many times with its starting longitude iterated by 500 m each time.This method allows high density coverage whilst maintaining the correct satellite perspective for realistic averaging.The footprints which lie within the boundary of the limits given in Table 1 are flagged and used to create the retrieval locations.
Each of the orbital positions given by the orbit simulator defines the location of the on-line footprint with the offline footprint being defined as the next vertically adjacent MODIS data value.
The error in the CO 2 retrieval is calculated by simulating a perfect (a-priori) retrieval (V = 0) and comparing this to a total column retrieval calculated using modified off-line surface reflectance's (equations 14, 15 and 16).
A bias is expected to be present in the simulated CO 2 retrieval error owing to the TC-DIAL equation being nonlinear and containing uncertainties in S on and S of f .Amediek et al. (2009) avoided the presence of a retrieval bias in their simulated retrievals by separating the measured differential optical depth into two terms (measurement and surface offset), and applying a power law expansion to convert the natural log of the surface ratio into a relative difference, therefore removing the possibility of retrieval error bias via the surface contribution term (equations 23 and 24).
∆τ gas is the atmospheric differential optical depth, S of f and S on are the measured off and on-line pulse energies, and ρ of f and ρ on are the off and on-line surface reflectance's (Amediek et al., 2009).
To provide continuity, the approach used by Amediek et al. (2009) was investigated for application in this study.It was found that the power law approximation was incompatible with the a-priori error approach adopted here.The information provided by the considered TC-DIAL instrument (S on and S of f ) does not allow discrimination of the reflectance variability contribution from the atmospheric contribution.It is therefore seen to be inappropriate in this instance to separate out the optical depth contributors so as to apply a power expansion to negate the presence of a retrieval bias.Furthermore, it is likely the inaccuracy of the power expansion will lead to some level of error in the optical depth and therefore the retrieval, especially in cases where the surface reflectance variability is high.Therefore, the Amediek et al. (2009) method was not used as a means of avoiding the retrieval error bias in this study.

Results
Each TC-DIAL measurement consists of multiple soundings integrated over a distance on the Earth's surface.A number of possible integration methods may be applied.The approach chosen is to calculate the arithmatic mean of the retrievals from the individual atmospheric soundings.The retrieval error distributions produced using this averaging process on the data obtained from each regional simulation are defined by their mean and standard deviations in Table 3 The error statistics in Table 3 indicate a significant regional variability in both the mean and standard deviations of the retrieval errors, most notably between America and Europe.The maximum regional difference in the mean and standard deviation is approximately 0.3 ppm and 0.1 ppm respectively.These large variations are indicative of the differences in the size and layout of agriculture, especially between Europe and America.
The retrieval averaging approach allows the highest level of information to be applied to the retrievals as atmospheric profiles are defined for each sounding.The addition of an error to the non-linear TC-DIAL equation introduces a bias which is relatively inconsequential for individual measurements (Figure 4) but significantly more important after the application of the averaging procedure (figure 5).
Fig. 4. Single sounding retrieval error prior to averaging for all simulated regions combined using the interpolation scaling method.Mean = 0.092 ppm, sd = 2.4 ppm Fig. 5. 50 km integrated single sounding retrieval errors for all simulated regions combined using the interpolation scaling method.Mean = 0.247 ppm, sd = 0.155 ppm An alternative process of integration which does not produce a bias is the averaging of the received laser intensities, allowing only one retrieval to be carried out per integration distance.The consequences of this method of averaging are the manifestation of extra uncertainties, such as the errors associated with the numerical weighting of surfaces with higher albedos over those with lower albedos, and the errors associated with the assumption that a single atmospheric profile is able to represent entire 50 km lengths of atmosphere.
Figure 6 is plotted to show an example of the application of the intensity averaging technique.The uncertainties in the retrievals that would result from such a procedure results in this method not being used here to derive error distribution statistics.
Fig. 6.50 km integrated retrieval error using received intensity averaging approach and interpolation scaling method.Produced as an example of having no bias.Mean = 0.002 ppm, sd = 0.093 ppm

Conclusions
A retrieval bias in the integrated measurements has been identified with an average magnitude of 0.25 ppm and a regional variability ranging from approximately 0.1 to 0.4 ppm across the scenes considered.Removal of the bias has not been possible in this study owing to the retrieval method of averaging and the a-priori method by which the errors were determined.The results obtained are representative of those that would be observed by a spaceborne TC-DIAL system should it's retrieval adopt the same averaging scheme as the one used here.
The results derived are constrained to specific regions of agriculture owing to the requirement that the surfaces considered are reasonably flat and shadow free.The importance of agriculture in the carbon cycle and its high surface reflectance variability has led to it being focused upon in this study, however the application of the methods described here are not necessarily limited to this biome alone.
A direct comparison with the results from Amediek et al. ( 2009) is relatively difficult owing to the difference in the range of surface types considered and the presence of a bias in this paper, however, the magnitude of errors from both papers are indeed similar, with Amediek et al. (2009) quoting an overall RMS error of 0.22 ppm.
It has therefore been demonstrated that spaceborne BRDF data may be downscaled using higher resolution surface radiance data to simulate the reflectance variability error in a TC-DIAL retrieval.Utilizing an a-priori approach, the actual error observed in a TC-DIAL retrieval has been calculated and a regionally varying retrieval bias has been shown to be present.Careful consideration for this regional variability would be required if such TC-DIAL measurements were to be assimilated into a carbon flux model.

Fig. 1 .
Fig. 1.Diagrammatic representation of a satellite based DIAL measurement showing four soundings within an integration interval.The two transmissions which make up a single measurement are closely located but do not completely overlap

Fig. 2 .
Fig. 2. Surface reflectance map over Europe generated by the Leicester LiDAR Model at 500m resolution

Fig. 3 .
Fig. 3. Landsat 7 reflectance data for two vertically adjacent MPAs with a single TC-DIAL viewing geometry LAT highlighted in colour

Table 1 .
500 m to 30 m and 10 m spatial resolution scaling factors for 5 agricultural scenes.

Table 2 .
List of model parameters used in the LLM for this study

Table 3 .
. Surface reflectance uncertainty distribution statistics for three scaling methods in units of ppm