Non-Lambertian Corrected Albedo and Vegetation Index for Estimating Land Evapotranspiration in a Heterogeneous Semi-Arid Landscape

Abstract

The application of energy balance algorithms to remotely sensed imagery often fails to account for surface roughness variation with diverse land cover, resulting in poor resolution of evapotranspiration (ET) variations. Furthermore, the assumption of a horizontally homogeneous Lambertian surface reflecting energy equally in all directions affects the calculations of albedo and vegetation index. The primary objective of this study is to improve the accuracy of the estimation and discrimination of ET among different land cover types in Southern New Mexico from ASTER datasets, by formulating the spatial variation of non-Lambertian reflectance using a wavelength-dependent Minnaert function.

1. Introduction

Evapotranspiration (ET) represents the water loss from the land surface to the atmosphere, as evaporation from soil or plant surfaces and transpiration from plants. Since water is a limiting resource in many ecosystems, it is important to estimate ET over the landscape accurately for planning and management of water resources, especially in arid or semi-arid regions such as Southern New Mexico. However, modeling the spatial variation of ET over heterogeneous landscapes is a difficult task for meteorologists, hydrologists, and agronomists [1].

ET is calculated from energy balance models, such as the Surface Energy Balance Algorithm for Land (SEBAL) [2], which is used in this study. It is based on the energy balance equation: Rn = H + G + λET, where Rn = net surface radiation including both short and long wave radiations (W/m²); H = sensible heat flux (W/m²); G = soil heat flux (W/m²); λET = latent heat flux (W/m²). These algorithms are based on data collected by satellites that provide information on surface temperature from TIR (thermal infrared) bands, on vegetation index (VI) from VNIR (visible near-infrared), and on albedo from VNIR and SWIR (shortwave infrared). They require little ground-based meteorological data (mainly wind speed, air temperature, and air humidity). Then, from the above input data, the values of net radiation, vegetation index, albedo, roughness length, soil heat flux, and sensible heat flux are calculated. Final ET is then derived as the residual term of the energy balance equation. For energy balance algorithm description, see [3,4,5].

The topographic effect on the radiance detected by satellite sensors is well known. Specifically, its influence on the solar illumination may affect the albedo, and cause misclassification of similar land cover areas that reflect the light differently. ASTER surface radiance and reflectance products are normalized by a simple cosine factor correction of the incident irradiance, in which the angle between the surface and the sun is computed by using a digital elevation model. According to Lambert’s cosine law [6], the amount of irradiance reaching an inclined pixel is proportional to the cosine of the incidence angle:

L_H = L_T (cos_sz/cos_i)

(1)

where L_H = radiance on a horizontal surface, L_T = observed radiance on a sloped surface, sz = solar zenith angle, i = sun incidence angle on the terrain. The solar illumination, cos_i, is calculated as:

cos_i = cos_sz cos_ts + sin_sz sin_ts cos_(sa−ta)

(2)

where ts = terrain slope, sa = solar azimuth, ta = terrain aspect. This topographic normalization is commonly used in energy balance models. It is based, for simplicity and convenience, on the unrealistic assumption that the earth is flat with a horizontally homogeneous Lambertian reflectance surface [7], which reflects energy equally in all directions. Also, it corrects only for differences in illumination caused by the orientation of the surface [8]. The bidirectional reflectance distribution function (BRDF) of real soil or vegetation deviates from Lambertian, and in a manner that varies with wavelength. In this way, images tend to be over-corrected, with slopes facing away from the sun appearing brighter than sun-facing slopes due to diffuse sunlight being relatively more influential on the shady slope [8,9]. Considering that most objects on earth show a non-Lambertian behavior, this model should be viewed with some skepticism for areas with slopes on the scale of several pixels.

VI and albedo (total short wave broadband reflectivity), which is simplified as ‘hemispherical’ under the Lambertian assumption, derive from arithmetic operations computed on surface reflectance bands of satellite imagery. Albedo is used to calculate the net radiation flux, and is derived from reflectance bands by commonly applying the regression model of Liang [10]:

α = 1.484α₁ + 0.335α₃ − 0.324α₅ + 0.551α₆ + 0.305α₈ − 0.367α₉ − 0.0015

(3)

where α_i is the spectral albedo for each ASTER band. The Modified Soil Adjusted Vegetation Index (MSAVI) is used to calculate the soil heat flux and the surface roughness, and derives from near-infrared and red bands: MSAVI = {(2NIR + 1) − [(2 NIR + 1)² − 8(NIR − red)]^1/2} / 2) [11]. Since albedo and VI are among the key factors affecting the energy balance calculation, errors in their estimation, will affect the calculation of ET [12]. Hence, in this study a non-Lambertian reflectance correction with the Minnaert function (Equation 4) [13] is applied to the ASTER surface reflectance bands to account for anisotropic surfaces among the different land cover types in the landscape.

L_H = L_T (cos_sz/cos_i)^K_α

(4)

where K is the Minnaert constant, and α = spectral band. The value of K depends on the optical properties of the ground cover, wavelength, and topography. It describes the roughness of the terrain and weakens the power of the Lambertian topographic correction. It has the advantage of being applied to all bands, and is obtained (Equation 5) by logarithmically linearizing Equation 4, and then estimating the slope of the regression between the radiometrically distorted bands and a shaded terrain model [14]:

log L_T = log L_H + K log (cos_i/cos_sz)

(5)

The objective of this study is to improve the accuracy of remotely-sensed ET calculations across different land cover types by analyzing the spatial variation of anisotropic reflectance and roughness using a non-Lambertian topographic correction on a pixel-by-pixel basis and processing each band separately. Then, the corrected albedo and VI will be applied to SEBAL for final ET calculation.

2. Methods

2.1. Study Area

The study area Jornada Experimental Range (JER) (~103,000 ha) is located in the Northern Chihuahuan Desert in Southern New Mexico, between the flanks of Mount Summerford on the southwest and the piedmont of the San Andres Mountains on the east, with an elevation range between 1,210–1,400 m (Figure 1). The climate is semi-arid, characterized by abundant sunshine, low relative humidity, and a total annual rainfall of approximately 24 cm, of which nearly 66% fall from May through October, primarily as intense, convective thundershowers [15]. The mean monthly maximum temperature is highest in June (36 °C). This semidesert area is a mosaic of shrub communities interspersed with grasslands comprising of five main habitats: creosotebush scrub (Larrea tridentata); tarbush shrubland (Fluorensia cernua); mesquite duneland including mesquite (Prosopis glandulosa, soil accumulations at base of plants 20 cm or less), mesquite dunes (20 cm to 3 m in height), and mesquite sandhills (greater than 3 m in height, sometimes 6 m or more); upland grassland dominated by blackgrama (Bouteloua eriopoda); and playa grassland (periodically flooded areas that receive drainage waters from upslope areas) dominated by tobosa (Pleauraphis mutica), dropseed (Sporobolus spp., threeawn (Aristida sp.), broom snakeweed (Gutierrezia sarothrae), and burrograss (Scleropogon brevifolius) [16].

Figure 1. ASTER surface reflectance VNIR band 3N over Las Cruces area on 08 June 2005.

Figure 1. ASTER surface reflectance VNIR band 3N over Las Cruces area on 08 June 2005.

2.2. Satellite, Terrain, and Meteorological Data

Source data encompass Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) level 2-B images of surface reflectance and surface kinetic temperature. The surface reflectance products consist of three visible near-infrared (VNIR) bands (0.52–0.86 μm) at 15 m resolution, and six short-wave infrared (SWIR) bands (1.60–2.43 μm) at 30 m resolution. The surface temperature products include five thermal infrared (TIR) bands (8–12 μm) at 90 m spatial resolution. All data are atmospherically corrected [17]. Acquisition date and time for the scene analyzed here was 8 June 2005, 17:56:02 (GMT) under cloudless conditions. A Digital Elevation Model (DEM; 7.5-min Quad 30-meter spatial resolution) was used for imagery pre-processing, derivation of slope and aspect map layers, and surface temperature calibration of dry adiabatic rate for choice of the hot and cold anchor pixels required for the SEBAL algorithm. A land cover shapefile of the study area containing information on polygon boundaries, plant species, and species dominance for each land cover type was used for reflectance correction. Solar zenith and azimuth angles were 80.5 and 171.6 (E of N) degrees respectively at time of satellite overpass (from Astronomical Application, Dept. US Naval Observatory, Washington, DC, USA). The solar radiation level, mean air temperature, mean wind speed, relative humidity, and vapor pressure for the time of satellite overpass were obtained from local climatological stations, and were used to calculate the near-surface temperature difference (surface temperature vs. air temperature). The last rainfall prior to the scene acquisition took place between 26 and 28 May with a total of ~19 mm averaged among the Long Term Ecological Research weather stations and the supplementary network of 11-rain gauges located along the study area elevation gradient. Ground measurements of latent heat flux at the time of satellite overpass were obtained from two eddy-covariance (EC) towers located on the upland grassland and mesquite dunes sites at the JER. The latent and sensible heat fluxes were corrected using the Webb-Pearman-Leuning equation [18] to account for density and heat effects. The distribution of EC ET values for the two sites is shown in Figure 2.

Figure 2. Eddy-Covariance ET flux data at the JER in June 2005.

Figure 2. Eddy-Covariance ET flux data at the JER in June 2005.

2.3. Technique and Data Layers Pre-Processing

ASTER imagery was integrated into a geographic information system (GIS). Data conversion, geoprocessing, and map algebra functions for spatial analysis computations of ET algorithms were performed in ArcGIS 9.3 software. ASTER images, DEM, and land cover coordinates were projected to the Universal Transverse Mercator (UTM) system, zone 13N, WGS 84 datum. The ASTER TIR, VNIR, and SWIR images were first orthorectified using the DEM, geodetic, and elevation control data to correct for positional accuracy and relief displacement. MSAVI was computed from ASTER reflectance bands 2 and 3N.

2.4. Spatial Computation of Minnaert Coefficients and Application of Reflectance Datasets to SEBAL Algorithms

A spatial analysis of the non-Lambertian behavior for each land cover type was performed on a pixel-by-pixel basis for the reflectance channels 1, 2, 3N, 5, 6, 8, and 9. Pixels of each raster reflectance band that corresponded to the areas defined by each vegetation class were extracted for a total of 91 raster datasets (7 bands × 13 land cover types). To calculate the Minnaert K coefficient for each of the extracted raster datasets, a backward reflectance correction calculation was computed using map algebra language. From the original ASTER reflectance cosine corrected values (Equation 1), L_T is calculated as a residual for each single pixel (Equation 4). A cos_i raster dataset was created from Equation (2), and a log of L_T (Equation 5) raster dataset was obtained for each vegetation type and each band. From Equation (5), the K coefficient is then calculated as the slope of the linear regression relating log L_T and log (cos_i/cos_sz). Regressions for each vegetation type were calculated using the ‘ArcView Grid and Theme regression’ extension [19]. The K values obtained were applied to Equation (4) for each land cover type and each band to obtain corrected reflectance datasets, which were then mosaicked into one final Non-Lambertian reflectance raster. Two surface albedo datasets were derived from the original Lambertian reflectance and the corrected non-Lambertian reflectance bands. Likewise, two sets of MSAVI were calculated. Albedo and MSAVI datasets were then resampled to 90 m resolution to match the ASTER surface temperature resolution, and were applied to the SEBAL energy balance algorithm. ET was calculated on a pixel-by-pixel basis at 90 m spatial resolution. Two final ET maps were generated for comparison, one using the Lambertian reflectance and one using the non-Lambertian reflectance.

To validate the ET models estimates, values were compared to ground fluxes data obtained from the eddy-covariance towers. A flux footprint model [20] was developed to locate the dominant ET source pixels. To assess the statistical power of discrimination of ET values between land cover types for the corrected and non-corrected models, a pairwise comparison of mean ET values was performed using the Scheffe’s test [21].

3. Results and Discussion

3.1. Minnaert-Corrected Reflectance

Each of the 13 land covers for each band was assigned a calculated Minnaert K coefficient (Figure 3). K varies between 0 and 1 for most surfaces but can be >1 when there is a dominance of the specular radiance component [22]; when K = 1, the surface is Lambertian. Playa grasses showed the lowest K values (0.64–1.05, mean = 0.91); while the upland grass blackgrama showed higher values (0.98–1.06, mean 1.02). Shrubs (creosote, tarbush, and ‘other shrubs’) and mesquite sandhills showed the highest values (0.994–1.07, mean = 1.03; and 1.01–1.04, mean = 1.02, respectively). Mesquite dunes showed values close to 1 (Lambertian surface; mean = 1.00). The values of the K factor vary across bands for upland and playa grasslands, and are similar for duneland and shrubs. Overall, they are higher for band 3 (0.76–0.86 µ; mean K = 0.98) and lower for band 1 (0.52–0.60 µ; mean K = 0.95). With K occurring as an exponent in the Equation 4, even small deviations of its value from 1 result in a considerable change of the equation’s solution. Then, the K values close to 1 obtained for the mesquite dunes still produce a change in the equation which propagates in the final calculation of ET, although this change was the smallest among the land cover types.

Figure 3. Minnaert regressions for each land cover class and each ASTER reflectance band at the JER.

Figure 3. Minnaert regressions for each land cover class and each ASTER reflectance band at the JER.

Compared to the original reflectance image, the Minnaert-corrected reflectance (Figure 4) and derived albedo resulted in higher variability of values among land cover types.

Figure 4. Comparison of reflectance (a) before and (b) after correction.

Figure 4. Comparison of reflectance (a) before and (b) after correction.

The playa grassland (illustrated in Figure 7) shows a distinctly lower radiance than the other land cover types due to the high-density of its canopy cover. The expected high reflectivity of sparse vegetated areas due to their sandy soil component was indeed found on tarbush shrubland, creosote scrubland, and duneland, especially mesquite sandhills. In contrast, the original reflectance showed similar values for playa grasses and sparse vegetated areas.

The effect of the Minnaert-corrected reflectance on the vegetation index is shown in Figure 5. Higher MSAVI values on grassland and lower values on tarbush, creosote, and sandhills better characterize the highly heterogeneous environment in vegetative cover and species composition of the JER.

Figure 5. Comparison of MSAVI (a) before and (b) after reflectance correction.

Figure 5. Comparison of MSAVI (a) before and (b) after reflectance correction.

3.2. Assessment of Spatial Distribution of ET

Derived ET maps from the corrected and non-corrected SEBAL are displayed in Figure 6.

Figure 6. Final (a) non-reflectance corrected and (b) reflectance corrected ET (mm/day) maps, which show a pattern with (c) surface temperature and (d) land cover type, respectively.

Figure 6. Final (a) non-reflectance corrected and (b) reflectance corrected ET (mm/day) maps, which show a pattern with (c) surface temperature and (d) land cover type, respectively.

As shown in Figure 6, the original SEBAL image (a) reveals a spatial distribution of ET that corresponds with the surface temperature distribution (c) across the area, as discussed below. The modification of the reflectance bands has generated a more heterogeneous spatial distribution of ET (b) that strongly matches with the different land cover types across the entire heterogeneous vegetated area (d). In the modified model, the lowest values were obtained over the mesquite sandhills and mesquite dunes, where the vegetation density is low, and higher values were obtained over fully vegetated surface of the playa grasslands, where ground depressions allow water drainage from upland slope areas and consequent water accumulation. The surface contours map is generated through a triangular irregular network or TIN (Figure 7). Over the sparsely vegetated surfaces of creosote scrubland and tarbush shrubland, intermediate values of ET were computed. Similar values were obtained on the upland blackgrama grassland, probably due to the water flowing to the downslope Eastern areas. This water accumulation supports both growth and transpiration of the grasses. The non-corrected SEBAL shows more homogeneous values across the study site. The strong correspondence of ET values with the surface temperature reveals a lower capability of modeling the surface resistance to heat transfer induced by the different vegetation covers. Aerodynamic resistance is underestimated over sparsely vegetated surfaces.

Figure 7. TIN from contours. Playa grassland is located at the lowest elevation of the area.

Figure 7. TIN from contours. Playa grassland is located at the lowest elevation of the area.

The scatterplot (Figure 8) of single ET pixel values shows poor agreement between the corrected versus the non-corrected SEBAL, with a high and low range of values, respectively. As shown in Figure 9, the non-corrected SEBAL produced overestimation of ET values over dunelands and shrublands, and underestimation on the playa grassland.

Figure 8. Scatterplot of ET single pixels values (mm/day) of corrected SEBAL (y) versus non-corrected SEBAL (x).

Figure 8. Scatterplot of ET single pixels values (mm/day) of corrected SEBAL (y) versus non-corrected SEBAL (x).

Figure 9. Scatterplots of ET single pixels values (mm/day) for (a) playa grassland, (b) creosotebush scrub, (c) tarbush shrubland, and (d,e) mesquite duneland for corrected model (y) and non-corrected (x).

Figure 9. Scatterplots of ET single pixels values (mm/day) for (a) playa grassland, (b) creosotebush scrub, (c) tarbush shrubland, and (d,e) mesquite duneland for corrected model (y) and non-corrected (x).

The power of both models to discriminate the spatial distribution of ET among the different land cover types was assessed using the Sheffe’s Test, as shown in Figure 10. The mean ET values were significantly different between 76 out 78 pairwise land covers comparisons (α = 0.05) in the corrected SEBAL (only differences between Aristida and tobosa, and between creosote and tarbush were not significant). In the non-corrected SEBAL, only seven pairwise comparisons were significantly different: blackgrama from tarbush, burrograss, and ‘other grasses’, and Aristida or burrograss from Sporobolus and mesquite dunes.

Comparison between direct ground measurements of ET by Eddy-Covariance towers (grass site ET = 38.98 W/m²; mesquite dunes site ET= 8.15 W/m²) and corrected SEBAL shows a good agreement at both sites. The ET value from the corrected-SEBAL at the grass site was 47.7 W/m², and at the mesquite dunes site was 11 W/m². The non-reflectance corrected SEBAL showed ET values of 99.7 W/m² at the grass site and 75 W/m² at the dunes site. At the mesquite dunes site, the EC ET fluxes are substantially lower relative to the grass site (Figure 2), clearly reflecting the real surface conditions of these two sites. In general, the comparison of large scale remote sensing data with ground-truth point measurements representing tiny areas is a known limitation in remote sensing studies including energy balance of ET. However, in this study similar SEBAL ET values to the EC ET values at the grass site were obtained in the other equal land cover areas. Likewise, similar SEBAL ET values to the EC ET values at the mesquite dunes site were obtained in the other mesquite dunes areas throughout the study site. Contrarily, the non-corrected model showed an overall similarity in ET values among the different land cover types.

Figure 10. Mean ET (mm/day) of each land cover type.

Figure 10. Mean ET (mm/day) of each land cover type.

4. Conclusions

This study demonstrates a new approach for improving the energy balance estimates of evapotranspiration using SEBAL, and its regional variability on a heterogeneous landscape. Since most vegetation canopies are non-Lambertian reflectors, a Minnaert correction of the original ASTER reflectance was performed to analyze the effect of the anisotropic behavior of the different land cover types on albedo and VI, and consequent effect on ET. The non-Lambertian response of the land cover, especially playa grassland, mesquite sandhills, creosotebush, tarbush scrubland, and other shrubs, generated a higher range of albedo and MSAVI values among the different land cover types in the corrected than in the non-corrected images. Considering that the Jornada Experimental Range is highly heterogeneous in vegetative cover and species composition, the reflectance-corrected MSAVI image better represents this pattern showing higher values for the highly vegetated grassland than the semi-arid shrubland and duneland. A comparison between the spatial distribution of ET in the corrected and non-corrected models confirms a higher heterogeneity of values of the former across the landscape with higher, intermediate, and lower ET values among grasslands, schrublands, and dunelands respectively, while the latter resulted in more homogeneous ET values among land covers and in overestimation of ET over dunelands and underestimation over grasslands. Also, the corrected model showed stronger agreement with the Eddy-Covariance measurements than the non-modified model. The different distribution pattern of the non-corrected model seems fundamentaly to reflect the surface temperature characteristics, and to inadvertently underestimate the aerodynamic resistance to heat transfer induced by the different roughness of the single vegetation types. Hence, the non-Lambertian topographic correction is a key factor in using reflectance for ET estimation in heterogeneous landscape.

In this study, a Minnaert coefficient was calculated for each land cover and each band from pixel-based slope and aspect raster datasets, without specifically evaluating the effect of specific slopes and aspects on the land covers radiance. Hence, future work should address a non-Lambertian coefficients computation based on a combination of land covers, slopes, and aspects.