Disaggregation of regional-scale (103 m) flux estimates to micrometeorological scales (101–102 m) facilitates direct comparison between land surface models and ground-based observations. Inversely, it also provides a means for upscaling flux-tower information into a regional context. The utility of the Atmosphere–Land Exchange Inverse (ALEXI) model and associated disaggregation technique (DisALEXI) in effecting regional to local downscaling is demonstrated in an application to thermal imagery collected with the Geostationary Operational Environmental Satellite (GOES) (5-km resolution) and Landsat (60-m resolution) over the state of Oklahoma on 4 days during 2000–01. A related algorithm (DisTrad) sharpens thermal imagery to resolutions associated with visible–near-infrared bands (30 m on Landsat), extending the range in scales achievable through disaggregation. The accuracy and utility of this combined multiscale modeling system is evaluated quantitatively in comparison with measurements made with flux towers in the Oklahoma Mesonet and qualitatively in terms of enhanced information content that emerges at high resolution where flux patterns can be identified with recognizable surface phenomena.
Disaggregated flux fields at 30-m resolution were reaggregated over an area approximating the tower flux footprint and agreed with observed fluxes to within 10%. In contrast, 5-km flux predictions from ALEXI showed a higher relative error of 17% because of the gross mismatch in scale between model and measurement, highlighting the efficacy of disaggregation as a means for validating regional-scale flux predictions over heterogeneous landscapes. Sharpening the thermal inputs to DisALEXI with DisTrad did not improve agreement with observations in comparison with a simple bilinear interpolation technique because the sharpening interval associated with Landsat (60–30 m) was much smaller than the dominant scale of heterogeneity (200–500 m) in the scenes studied. Greater benefit is expected in application to Moderate Resolution Imaging Spectroradiometer (MODIS) data, where the potential sharpening interval (1 km to 250 m) brackets the typical agricultural field scale. Thermal sharpening did, however, significantly improve output in terms of visual information content and model convergence rate.
The set of instruments deployed on the current fleet of satellite remote sensing platforms provides unique possibilities for monitoring the biotic and hydrologic condition of terrestrial ecosystems and their response to environmental changes. In particular, these instruments allow mapping of evapotranspiration (ET) and moisture stress at a wide range in spatial resolution, covering areas from local to regional to global scales. The breadth of resolution afforded by these combined instruments is crucial to engendering a broad understanding of land surface response. High-resolution (101 m) data can help to pinpoint the cause of remotely sensed phenomena, while lower resolution (104 m) imagery are better suited for assessing impacts over large areas.
Norman et al. (2003) presented a nested flux-disaggregation technique that has demonstrated utility in integrating multiscale remote sensing data to effect upscaling and downscaling in surface fluxes (see Fig. 1). Disaggregation provides a potential means for bridging extant gaps in spatial scale that have traditionally posed challenges to surface-flux-monitoring programs: physical gaps encountered in upscaling ground observations to landscape and regional scales (or, equivalently, downscaling regional model predictions to the ground-observation footprint scale) and logistical gaps in coverage provided by current satellites, which have either high-temporal/low-spatial resolution or low-temporal/high-spatial resolution.
New-generation flux-monitoring campaigns are recognizing the need to establish connectivity between intensive surface observations, typically collected at a point or a distributed set of points, by placing them within a regional context using modeling supported by remote sensing (Running et al. 1999). The North American Carbon Program, for example, is adopting a multitiered sampling scheme (Wofsy et al. 2002) in which general surface biophysical and land-use properties will be monitored over continental scales using remote sensing, with more intensive sampling of carbon stores and fluxes occurring on the ground at sites selected to represent the endemic range in spatial variability. Large-scale carbon flux networks, such as AmeriFlux (Baldocchi et al. 2001) and EuroFlux (Valentini 2003), will require robust methodologies for upscaling and integrating observations made at individual towers to be able to draw regional inferences regarding terrestrial carbon cycles (Fan et al. 1998).
From a modeling standpoint, the regional-scale soil–plant–atmosphere models required to support these upscaling efforts are notoriously difficult to validate because of the gross scale mismatch between the model grid cell (typically 10–100 km) and the surface footprint sampled by a ground-based micrometeorological tower (0.1–1 km). The problem of validation is exacerbated over heterogeneous landscapes, such as those depicted in Fig. 1, where great care must be taken to select representative ground sampling sites. Improved estimates of mean surface flux conditions can be obtained by averaging fluxes sampled at multiple locations within a model grid cell (e.g., Doran et al. 1998; Gao et al. 1998) or along an aircraft transect (Desjardins et al. 1997; Oechel et al. 1998; Mahrt et al. 2001; Kustas et al. 2001b; Chen et al. 2003); however, such measurements typically require large investments of time and resources and therefore cannot be widely implemented. Alternatively, methods capable of disaggregating regional flux estimates down to the subkilometer scale provide a means for quantitative validation against observations made at individual towers.
Large-scale remote sensing projects are also confronted with the temporal/spatial coverage trade-off inherent in the current suite of earth-observing satellites. While low-spatial-resolution sensors like the Geostationary Operational Environmental Satellite (GOES) can provide good temporal coverage (15 min), high-spatial-resolution instruments, as on the Land Remote Sensing Satellite (Landsat), tend to have more limited temporal coverage (biweekly) and yield little information about diurnal variability in land surface behavior. In addition, the thermal infrared (TIR) bands, which can provide important information about surface moisture status, are typically instrumented at significantly lower spatial resolution than are the visible–near-infrared (vis–NIR) bands, which can be used to estimate vegetation-cover amount (see Table 1). For many applications, synergistic use of both high-temporal-resolution and high-spatial-resolution datasets in multiple wave bands can be advantageous.
Methods have been developed for combining high-resolution vis–NIR and TIR remote sensing data to deduce the surface energy balance, each with strengths and weaknesses in terms of routine application. Some rely on the availability of contemporaneous in situ measurements, primarily near-surface meteorological conditions such as air temperature, wind speed, and humidity, and are therefore difficult to implement operationally over large regions (Gardner et al. 1992; Choudhury et al. 1994; Moran et al. 1994, 1996). Others require that a sufficient range in vegetation cover and surface temperature (Gillies and Carlson 1995; Jiang and Islam 2001) or hydrologic (dry–wet) conditions (Bastiaansen et al. 1998a,b) be present within the scene so that end points in the temperature–vegetation-cover relationship can be defined.
The nested-scale flux-modeling system introduced by Norman et al. (2003) uses high-resolution thermal, visible, and near-infrared imagery to disaggregate 5-km flux estimates from the regional Atmosphere–Land Exchange Inverse (ALEXI; Anderson et al. 1997) model down to the scale of a surface flux-tower footprint (101–102 m) without the need for any locally made measurements or restrictions regarding in-scene variability. This procedure [Disaggregated ALEXI (DisALEXI)] was applied to remote sensing data collected by aircraft at 24-m resolution over a 4-day dry down during the 1997 Southern Great Plains field experiment (SGP97). Footprint-weighted latent heat flux estimates agreed with eddy covariance measurements made at four flux towers to within 12% and captured temporal effects of the dry-down event (Norman et al. 2003).
As an ancillary enhancement to the DisALEXI technique, Kustas et al. (2003) describe an algorithm that uses the physical relationship between surface temperature and vegetation cover to sharpen TIR remote sensing imagery to higher resolutions associated with vis–NIR bands (Table 1). While the interband resolution difference is only a factor of 2 for Landsat 7 (30 versus 60 m), the discrepancy for the Moderate Resolution Imaging Spectroradiometer (MODIS), which provides better temporal coverage, is significant (250 m versus 1 km). Townshend and Justice (1988) argue that spatial resolutions of 250–500 m are required to effectively monitor disturbances in surface energy balance due to land-use and land-cover changes; thus, methodologies for improving important remote sensing inputs to such resolutions will be of value. Kustas et al. (2003) applied this sharpening algorithm (called DisTrad) to aircraft-acquired thermal imagery from SGP97, artificially degraded in resolution to simulate various satellite-band combinations listed in Table 1. DisTrad reconstructed high-spatial-resolution structure found in the native-resolution thermal imagery with good accuracy.
In this paper, the DisALEXI and DisTrad techniques will be applied to satellite data collected with the Enhanced Thematic Mapper Plus (ETM+) instrument on Landsat 7. Four Landsat scenes from 2000–01 are examined, each containing eddy covariance towers associated with the Oklahoma Mesonet (Brock et al. 1995). In addition to using satellite (as opposed to aircraft) data, this experiment covers a wider range of climatic and phenological conditions than those considered by Norman et al. (2003) and Kustas et al. (2003), who focused on 4 consecutive days in a small area around El Reno, Oklahoma. Disaggregated fluxes are compared to ground-based eddy covariance measurements, and spatiotemporal patterns in the distributed flux predictions are evaluated for reasonability. This study also examines the benefits of the DisTrad thermal-sharpening algorithm in terms of elucidating small-scale structure in modeled surface flux distributions and improving model behavior. The paper begins by briefly reviewing the ALEXI, DisALEXI, and DisTrad models.
2. Model descriptions
The ALEXI/DisALEXI modeling framework, diagrammed schematically in Fig. 2, is intended for diverse, routine applications and therefore attempts to balance the competing demands of generality and simplicity. The models have been designed to accommodate varying surface conditions while remaining computationally inexpensive and requiring only a tractable array of surface parameters that can be estimated remotely. To improve robustness, they are built where possible around differential (as opposed to absolute) radiometric signals (Kustas et al. 2001a).
a. The two-source model
At the core of this modeling system is a two-source (plant and substrate), land surface representation coupling conditions inside the canopy to fluxes from the substrate (typically soil), plants, and atmosphere. The two-source model (TSM; Norman et al. 1995b; Kustas and Norman 1999a,b, 2000) partitions the composite directional radiometric temperature [TRAD(ϕ)] of a heterogeneous scene into soil and canopy contributions (Ts and Tc) given an estimate of fractional vegetation cover apparent at thermal view angle ϕ [f(ϕ)]:
Equation (1) is solved along with a system of surface energy balance equations:
where RN is net radiation; H, LE, and G are fluxes of sensible, latent, and soil heat conduction, respectively; and the subscripts “c” and “s” represent fluxes from the soil and canopy components of the scene. Net radiation is estimated from modeled or measured downwelling shortwave and longwave radiation (upwelling components are modeled by the TSM), while G is related to the net radiation just above the soil surface (RNs), following Norman et al. (2000). Given specification of upper boundary conditions in air temperature and transport resistances based on wind speed, vegetation cover, and surface roughness, Hc and Hs are computed assuming the series resistance network in Fig. 2. A modified Priestley–Taylor (1972) relationship provides an initial estimate of canopy evapotranspiration (LEc), and the soil evaporation rate (LEs) is computed as a residual. If LEs is negative during the day (condensation onto the soil), this is taken as a sign of moisture stress, and LEc is throttled back.
The two-source representation is a major improvement over previous single-layer thermal models that required site-specific adjustments to compensate for differences in aerodynamic coupling between the soil, canopy, and atmosphere (Kustas et al. 1989; Hall et al. 1992; Stewart et al. 1994; Kubota and Sugita 1994). The TSM also provides a means for accommodating the dependence of apparent surface temperature on view angle [Eq. (1)], caused by the variable obscuration of the underlying bare soil when a canopy is viewed off nadir (Vining and Blad 1992). This feature is critical to remote sensing applications, where a scene may be viewed from many different angles.
Whereas lower boundary conditions are supplied by thermal remote sensing data, the TSM requires specification of temperature above the canopy and is particularly sensitive to biases in this input with respect to the TIR reference (Zhan et al. 1996; Anderson et al. 1997; Kustas and Norman 1997). Just a 1-K error in the assumed surface-to-air temperature difference can translate into errors in predicted sensible heating of up to 100 W m−2 (Norman et al. 1995a). Because shelter-level atmospheric properties can be strongly coupled to local surface fluxes, these upper boundary conditions generally cannot be interpolated with adequate accuracy from synoptic weather network observations, with a typical spacing of 100 km. For regional-scale flux-mapping applications, the TSM has been coupled with a simple model of atmospheric boundary layer (ABL) development (McNaughton and Spriggs 1986) so that near-surface air temperatures are simulated and consistent with modeled fluxes.
The Atmosphere–Land Exchange Inverse model is a coupled TSM–ABL model (Anderson et al. 1997; Mecikalski et al. 1999). The lower boundary conditions for ALEXI are provided by TIR observations taken at two times during the morning, separated by about 4 h, from a geostationary platform such as GOES. The ABL model relates the modeled rise in air temperature above the canopy during this interval and the resulting growth of the ABL to the time-integrated influx of sensible heating from the surface. This results in a prediction of spatially distributed fluxes at the time of the second surface temperature observation, about an hour before local noon.
Input data required for regional application of the ALEXI model are listed in Table 2 and reviewed by Mecikalski et al. (1999). Use of time-differential TIR measurements reduces model sensitivity to errors in absolute temperature due to sensor calibration and atmospheric and surface emissivity corrections. Importantly, the air temperature in the surface layer is not defined as a boundary condition—it is evaluated by the model at the TSM–ABL interface and responds to feedback from both the surface fluxes and the atmospheric profile. Sensitivity studies with ALEXI have shown that the modeled air temperature adjusts to instrumental biases in the TIR data, tending to preserve the true surface-to-air temperature gradient (Anderson et al. 1997).
Using local-scale TIR measurements made with ground-based infrared thermometers, which sample an area on the order of 10 m × 10 m, ALEXI flux predictions agree to within 20% with tower measurements made over a variety of land-cover types including corn, soybean, tallgrass prairie, sparse desert shrubs, and rangeland (Anderson et al. 1997, 2003). In practice, however, ALEXI is more suitably applied to satellite-based thermal data acquired at the 5–10-km scale—the scale at which organized land surface behavior becomes effective in influencing mean atmospheric conditions and driving boundary layer growth. To compare predictions at this scale with ground-based flux measurements, the regional-scale model fluxes need to be spatially disaggregated.
The DisALEXI algorithm (Norman et al. 2003) is a two-step process (Fig. 2): first, ALEXI is executed at a resolution of 5 km to estimate the air temperature at the interface between the land surface and ABL submodels [∼50 m above ground level (AGL)]; then, the TSM is applied to high-resolution surface temperature and vegetation-cover data (from aircraft or satellite), holding the air temperature at 50 m constant at the ALEXI-derived value. To ensure the ALEXI-derived vertical temperature gradient is preserved in the disaggregation stage, the high-resolution radiometric temperature field is corrected for biases (due to differences in sensor bandwidth, calibration, and/or atmospheric corrections) to yield an average consistent with the 5-km temperature data used in ALEXI.
The resulting high-resolution flux predictions can then be reaggregated and compared directly with tower observations, assuming some weighting function describing the surface source footprint contributing to the flux detected at the height and location of the sensor. Current techniques for estimating flux footprint distributions fall into two major categories: analytical solutions to the diffusion equation and Lagrangian models, which numerically simulate air parcel trajectories from source to detector. Analytical models (e.g., Schuepp et al. 1990; Horst and Weil 1992; Schmid 1994) use various simplifying approximations regarding turbulence parameterization and source variability but are relatively easy to implement. The more rigorous Lagrangian techniques (e.g., Leclerc and Thurtell 1990; Finn et al. 1996) require detailed specifications of the turbulence field, which may not be available in many applications.
The DisALEXI approach, like other “tile” or “mosaic” land surface models, assumes that horizontal fluxes are small in comparison with vertical fluxes and that conditions at 50 m AGL [approximating the atmospheric “blending height” (Wieringa 1986; Mason 1988) are more or less uniform over scales of 5 km. Both assumptions may be violated in strongly heterogeneous landscapes, where advection and small-scale surface–atmospheric coupling can become significant. The impacts of these assumptions under varying states of heterogeneity, and possible operational correction techniques, are being investigated in large-eddy simulations with an embedded TSM land surface representation (Albertson et al. 2001; Kustas and Albertson 2003).
d. Thermal sharpening
The ALEXI/DisALEXI models use remote sensing information from both the TIR (surface temperature) and vis/NIR (vegetation cover) wave bands. Model output resolution is typically limited by the resolution of the thermal sensor (∼102–103 m), which is often lower than would be desirable for local analyses and effective land surface change detection (∼250 m; Townshend and Justice 1988). On current satellites, vis–NIR sensors operate at 1–6 times higher resolution than thermal sensors (Table 1) and more closely approach this optimal resolution. The DisTrad algorithm developed by Kustas et al. (2003) sharpens the resolution of thermal band data to that of vis/NIR bands, taking advantage of the functional relationship between surface temperature and vegetation indices (VIs), such as the Normalized Difference Vegetation Index (NDVI; Rouse et al. 1973), observed within localized areas.
The DisTrad technique is based on fitting a second-order polynomial between radiometric temperature (TRAD the dependent variable) and the VI, aggregated to the coarser thermal resolution through linear averaging (VIlow, the independent variable):
Here, the “hat” symbol indicates a temperature value predicted by the resulting regression equation. Coarse-scale pixels used in this regression are selected from the scene such that each shows small subpixel heterogeneity (in terms of its coefficient of variation in VI), yet collectively represent a wide dynamic range in VI (see Kustas et al. 2003). Pixels with standing water are not considered.
The least squares relationship developed at the coarse scale [Eq. (3)] is then applied to the full-resolution VI data to predict temperature at that finer scale:
is a residual correction ensuring the average temperature in each coarse-scale pixel area matches the unsharpened temperature, TRAD,low. The restoration of these residuals, ΔT̂RAD,low, to the sharpened image is critical because it preserves spatial patterns in surface temperature that are due to moisture (rather than cover) variability and therefore are not accounted for by Eq. (3). The result is that high-spatial-frequency structure is created within each coarse-scale pixel such that it is consistent with overall temperature versus VI behavior exhibited within the scene, yet with the constraint that reaggregation will recover the original low-resolution thermal image.
In this study, the DisALEXI and DisTrad algorithms were applied to Landsat imagery collected over the state of Oklahoma, which hosts an extensive network of weather and flux-monitoring stations known as the Oklahoma Mesonet (Brock et al. 1995; Shafer et al. 2000).
The mesonet became operational in 1994 and consists of 115 automated stations that measure a set of core parameters including air temperature and relative humidity at 1.5 m, wind speed and direction at 10 m, atmospheric pressure, incoming solar radiation, rainfall, and bare and vegetated soil temperatures. In 1999, as part of the Oklahoma Atmosphere Surface-layer Instrumentation System (OASIS) Project (Brotzge et al. 1999), 90 mesonet stations were augmented with instruments for measuring surface heat fluxes using a gradient-profile technique (Paulson 1970). For comparison, eddy covariance instrumentation was additionally installed at 10 of these OASIS sites, designated as “supersites.” The analyses here focus on seven supersites that were located within Landsat 7 scenes acquired on 4 exceptionally clear days during 2000–01 (Fig. 1; Table 3).
a. Model input
1) Meteorological inputs
Ancillary surface and atmospheric data required by the modeling system include an estimate of the wind speed field at 50 m AGL and an early-morning atmospheric profile of temperature and mixing ratio (to 5–8-km altitude) at each 5-km grid cell in the ALEXI modeling domain (see Table 2). These input fields are currently created with the analysis component of a mesoscale model (in initialization mode) using standard observations from the synoptic weather and radiosonde networks. The temperature profile is used in the ABL submodel component of ALEXI; temperature and mixing ratio profiles are used to atmospherically correct the remotely sensed TIR imagery. Note that model sensitivity to errors in these meteorological inputs is relatively small (Anderson et al. 1997); therefore, local measurements are not required.
2) Coarse-resolution remote sensing data
The coarse-resolution surface temperature observations used in ALEXI were obtained with the GOES-8 imager instrument within the 10.2–11.2-μm (band 4) window. Thermal imager data are available every 15 min at an average nadir-viewing angle of approximately 40° and a nominal spatial resolution of 5 km at the location of the Oklahoma study area. Atmospheric corrections were performed using the methods described by French et al. (2003), and a vegetation-cover-dependent correction for surface emissivity was applied as in Mecikalski et al. (1999).
Downwelling solar and longwave radiation were estimated at each pixel in the ALEXI grid from hourly GOES data at 20-km resolution (Diak et al. 1996, 2000). Using a biweekly composited NDVI product (Eidenshink 1992) generated with the Advanced Very High Resolution Radiometer (AVHRR), estimates of vegetation-cover amount were derived with the expression suggested by Choudhury et al. (1994). These cover estimates were used in conjunction with a land surface classification created for the conterminous United States (USGS 1995) to assign relevant surface parameters such as surface roughness, displacement height, and radiometric properties (see Mecikalski et al. 1999 for more details).
3) High-spatial-resolution remote sensing data
High-resolution surface temperature and vegetation-cover maps were derived from Landsat7 ETM+ thermal and vis/NIR imagery. From Landsat archives for the years 2000–01, four predominantly clear scenes were selected, each containing two OASIS supersites (Table 3). A 5 km × 5 km area coincident with the ALEXI grid cell containing the supersites was then carved from these Landsat images, resulting in eight subscenes composed of 167 × 167 pixels of dimension 30 m.
The TIR data were extracted from ETM+ band 6 imagery, which are acquired at 60-m resolution. Each thermal subscene was atmospherically corrected with the method of French et al. (2003) using temperature and mixing ratio profiles from the corresponding ALEXI 5-km grid cell, and cover-dependent emissivity corrections were applied. Solar and longwave radiation estimates were also extracted from the ALEXI input fields and were held constant across each subscene.
NDVI was computed using 30-m resolution imagery from ETM+ bands 3 (visible) and 4 (NIR) and was used to map vegetation-cover fraction following Choudhury et al. (1994). To assign surface roughness characteristics, a land classification of “perennial ground cover” (class 7 in Mecikalski et al. 1999) was used uniformly across each subscene. This class was chosen to best represent the grassy land cover typical at OASIS tower locations, but will cause surface roughness to be underestimated in forested areas, resulting in locally depressed sensible heating rates. In future studies, classification maps should ideally be generated from high-resolution multispectral data.
b. Validation data
Eddy covariance (EC) flux measurements were collected at the seven OASIS supersites listed in Table 3. The EC instrumentation at these sites is described in detail by Brotzge et al. (1999), Brotzge (2000), and Basara and Crawford (2002). Sensible and latent heat fluxes are estimated with eddy covariance analysis and standard corrections using data from a 3D sonic anemometer and Krypton hygrometer mounted on an instrument tower at 4.5 m AGL. Additionally, net radiation is monitored with a four-component net radiometer mounted at 2 m AGL, while the soil heat flux is computed with the combination approach of Tanner (1960) using measurements from heat flux plates installed at 5-cm depth near the base of the tower, with heat storage in the surface soil layer estimated with two platinum resistance detectors at 0- and 5-cm depths. At each site, the instrument set is enclosed within a fenced area to deter vandalism and damage by grazing animals.
Errors in flux comparisons due to the difference between modeled and measured flux reference heights (50 and 4.5 m AGL, respectively) are expected to be small under the conditions studied here. Under clear noontime skies, and with sensible heating rates greater than 100 W m−2, the boundary layer depths are likely on the order of 1 km. Studies conducted in similar conditions and landscapes in Oklahoma concerning the use of 2-m-AGL versus mixed-layer meteorological measurements in the TSM indicate that the resulting differences in flux predictions were small (Kustas et al. 1999). Furthermore, MacPherson et al. (1999) found that aircraft flux observations at 25–35 m during SGP97 agreed well with 2-m flux-tower data when segmented by land use.
As is often the case with eddy covariance datasets (Kizer and Elliott 1991; Fritschen et al. 1992; Stannard et al. 1994; Lloyd et al. 1997; Twine et al. 2000; Wilson et al. 2002; Brotzge and Crawford 2003), the sum of the turbulent sensible and latent heat fluxes, H and LE, measured at the supersites was consistently less than the available energy, RN − G. Closure in the surface energy balance represented in these measurements ranged from 95% to <40%, as quantified by the ratio of (H + LE)/(RN − G). To maintain comparability with the models (which enforce closure), the flux measurements were corrected by modifying H and LE such that, over hour averaging intervals, they summed to the available energy (RN − G) yet retained the observed Bowen ratio (Twine et al. 2000). Fluxes from the BESS supersite consistently showed closure of less than 40%; therefore, measurements of LE and H from this site are not used in any quantitative flux comparisons.
Local meteorological datasets were also assembled for each supersite, although these data are used for interpretation purposes only. Time traces of precipitation and saturation deficit measured at the supersites are shown in Fig. 3 for a period of 2 weeks prior to each study day. Observations of wind speed and direction are used in the footprint analyses presented below.
a. Flux disaggregation
For each day of Landsat imagery, the ALEXI model was executed at 5-km spatial resolution over the entire Oklahoma modeling domain using a GOES-derived surface temperature rise signal measured from an hour past local sunrise to the time of the Landsat overpass (tLS), typically 4 h later. Output included modeled fields of the four major flux components and air temperature at 50 m AGL, each field coincident in time with the Landsat imagery. Modeled air temperatures were extracted from each ALEXI grid cell containing a target OASIS supersite and used as an upper boundary for disaggregating the 5-km fluxes associated with that cell (Fig. 2).
The atmospherically corrected Landsat TIR data used in the disaggregation were normalized with respect to the GOES temperature data to preserve the surface-to-air temperature gradient predicted by the ALEXI model. To minimize effects of errors in registration between the coarse-scale and finescale imagery, this normalization was determined over a 35 km × 35 km area (7 × 7 ALEXI cells) surrounding each mesonet station, yielding an average correction of +0.9°C. The bias-corrected 60-m thermal data were then interpolated to a 30-m grid coincident with the Landsat vegetation-cover data. For comparison, three different interpolation schemes were employed: a bilinear interpolation and two variations of the thermal-sharpening technique, DisTrad, as described below.
DisALEXI was applied to the resulting 30-m temperature and vegetation-cover fields for each Landsat subscene. As with ALEXI, the DisALEXI output flux fields represent a snapshot of instantaneous surface conditions at the time tLS. Jackson et al. (1983) suggested a procedure for extrapolating daily integrated surface fluxes from single-time estimates by assuming that the surface evaporative fraction (EF; the ratio of latent heat to the available energy) remains approximately constant throughout the day. Daily ET can then be related directly to the diurnal available energy curve, which can be predicted using only remote sensing and standard weather data. This is often a good approximation but can neglect effects due to afternoon changes in atmospheric conditions and/or surface moisture/temperature stress (Sugita and Brutsaert 1991; Brutsaert and Sugita 1992; Hall et al. 1992; Kustas et al. 1993).
For comparison with mesonet EC measurements, 30-m disaggregated fields of net radiation, and sensible, latent, and soil heat flux were reaggregated using the one-dimensional analytical footprint model of Schuepp et al. (1990, 1992), applied as a weighted average along a transect (∼30 m wide) at the angle of the mean local wind direction (averaged over the interval tLS ± 0.5 h). This weighting function uses estimates of local surface roughness and friction velocity, which were computed as weighted averages from model output along the same transect. While this is a relatively simplistic representation of complex, time-dependent transport phenomena, it does capture some signature of the azimuthal variability in surface conditions.
b. Thermal sharpening
The two-source model embedded in ALEXI and DisALEXI has difficulty converging given incongruous surface temperature and vegetation-cover data. Simple techniques for interpolating thermal data to the higher resolution of the cover information can result in incompatible inputs in localized areas, especially along narrow discontinuities such as roads and riverbeds. This leads to fringes of nonconvergent grid cells; for example, where the high radiometric temperature of an unresolved black asphalt road spills over into the full green cover in the roadside ditch. In such circumstances, model iteration using Eq. (1) can produce diverging estimates of soil temperature and sensible heating rates. The DisTrad thermal sharpening algorithm has the potential of improving the performance of the disaggregation process by interpolating structure into the thermal maps that is functionally consistent with the high-resolution cover information. The alternative is settling for the lower native resolution of the thermal sensor.
The sharpening tests presented here used NDVI as the vegetation index forming the independent variable in Eq. (3). Two methods of treating the residual term ΔT̂RAD,low in Eq. (4) were examined. Again, these residuals are critical to maintaining fidelity with the observed thermal data at their coarse native scale. One method is to restore the residuals in an unadulterated form, as implemented by Kustas et al. (2003); this will be referred to as the “standard sharpening method” (SS). This assures perfect correspondence with the low-resolution image but can give a boxy structure to the resultant sharpened image. The restored 60-m residuals often mask high-resolution structure reconstructed through the sharpening process.
A close examination showed that this boxiness could be removed, with little damage to coarse-scale fidelity, if the residuals are smoothed before being added back into the sharpened image. An optimal balance between fidelity and visual information content was achieved when the residual map was convolved with a Gaussian kernel with full width half maximum equal to the raw thermal resolution (60 m in this case). This will be referred to as the “convolved sharpening method” (CS).
Figure 4 demonstrates the utility of these thermal-sharpening techniques as applied to the Landsat subscene containing the IDAB mesonet site. Here, the raw 60-m Landsat temperature data have been interpolated to a 30-m grid using the standard and convolved-residual sharpening algorithms. For comparison, a simple bilinear interpolation (BL) from 60- to 30-m resolution is also shown. The convolved sharpening method is able to pull out significant detail in comparison with the other methods, aiding visual interpretation of features in the image. On the 60-m scale, however, similarity is maintained between the four thermal images.
5. Results and discussion
Disaggregated fields of latent heating surrounding the OASIS supersites are shown in Fig. 5. These maps demonstrate the influence of various thermal interpolation techniques (bilinear interpolation and standard and convolved sharpening) on model flux distributions. In the following subsections, the accuracy and utility of the disaggregation process and associated thermal-sharpening techniques are evaluated quantitatively in comparison with ground-based tower-flux data and qualitatively in terms of enhanced information content that emerges at high resolution where flux patterns can be identified with recognizable surface phenomena.
a. Model–tower-flux comparisons
Figure 6a compares instantaneous fluxes predicted by the ALEXI model at the 5-km scale with measurements made at the mesonet stations listed in Table 3. In comparison, Fig. 6b shows the improvement in agreement obtained when the 5-km fluxes are disaggregated down to the 30-m scale and integrated over the tower-flux footprint. Statistical measures of model performance in predicting the observed flux components RN, H, LE, and G are listed in Table 4, including the coefficient of efficiency (E) proposed by Nash and Sutcliffe (1970) as a performance metric preferable to the coefficient of determination (R2), which can indicate perfect agreement even in the presence of systematic model biases. The coefficient of efficiency is particularly useful in comparing the accuracy of predictions from multiple models with respect to a given set of observed data. With this index, the relative utility of various thermal-sharpening algorithms can be assessed in terms of DisALEXI flux prediction accuracy.
The benefit of using a disaggregation technique to validate kilometer-scale flux predictions made over heterogeneous landscapes is readily apparent in Fig. 6 and Table 4. While the overall bias in the ALEXI 5-km predictions is small in comparison with ground-based observations, the scatter is relatively large (18%) and will be strongly dependent on the exact location of the observing station within each ALEXI grid cell (see Fig. 5). The standard deviation in disaggregated latent heating estimates over a given 5 km × 5 km scene ranges from 40 to 80 W m−2; comparable to the root-mean-square difference (rmsd) of 90 W m−2 between ALEXI predictions and measured ET fluxes in Fig. 6a. To directly validate ALEXI predictions with tower data, multiple stations would need to be deployed within each cell, carefully positioned to capture a representative range in surface conditions. Disaggregation reduces the relative model error by a factor of 2 (to 9%) and yields a higher coefficient of efficiency in comparison with tower-based flux measurements.
The statistics in Table 4 indicate, however, that the specific technique used to interpolate the thermal Landsat data from 60 to 30 m did not significantly affect model accuracy. A simple bilinear interpolation performed as well as the standard and convolved thermal-sharpening algorithms over this range in resolution. This is consistent with the findings of Kustas et al. (2003), who note that flux predictions from DisALEXI were greatly improved using data sharpened to the 200-m scale, but little was gained by further sharpening to the 30-m scale. In that experiment, the dominant scale of temperature variations was defined by the typical size of the fields within the study scene (∼200 m); thus, sharpening to finer scales did not reap significant benefit. Similar scales of heterogeneity are evident in the scenes studied here.
The supersite showing the greatest discrepancy with model predictions is ALV2, where instantaneous latent heating is underestimated by 20% (∼90 W m−2). This site lies in pasture due north of a fenced field with low vegetation cover on the day in question (see Fig. 5a). High surface temperatures suggest the soil surface in this field had dried significantly since the heavy rainfall 3 days prior (Fig. 3). At the time of observation, winds at the site were from the south-southwest, causing the analytical footprint transect to intersect this field and lowering the footprint-integrated model latent heating. Proximity to the sharp discontinuity between pasture and field makes this site particularly sensitive to errors in image navigation and footprint analysis.
Consistent overestimation of G by both the ALEXI and DisALEXI models, evident in Fig. 6, may be due in part to subgrid-scale anomalies in vegetation cover coincident with the supersites. Most of the mesonet sites considered in this study are located in pasture, screened from animal intrusion by fencing. At some sites, the disparity in green cover fraction inside and outside of the site fencing is quite significant. The soil heat flux plates are buried within the protected area and therefore are susceptible to this type of small-scale cover bias, while the eddy covariance sensors integrate fluxes emanating from inside and outside the fencing. The upwelling radiation sensors on the tower-mounted net radiometer will similarly sample only a small area (∼10 m across) within the enclosure. In this study, however, the mean bias in modeled RN (8 W m−2) is small compared to that for G (33 W m−2). This is consistent with findings by Brotzge and Crawford (2003), who observed large differences (factor of 2) between diurnal soil heat curves measured at two flux stations spaced 100 m apart at the Foraker, Oklahoma, OASIS supersite (not studied here), but only slight differences in net radiation, suggesting that G may be more sensitive than RN to spatial variations in vegetation cover, soil type, and moisture content. Because midday net radiation is heavily driven by the shortwave component, particularly under clear-sky conditions, and because the albedo of the soil and vegetation is similar, this is not unreasonable.
Given the fact that the soil heat flux plates and eddy covariance flux systems are sampling different surface footprints, one might ask whether it is reasonable to use point-scale measurements of G to close the observed energy budget (Lloyd et al. 1997; Brotzge et al. 1999; Baldocchi 2003). A remote sensing model may in fact predict an average soil heat flux that is more representative of the source area contributing to the EC fluxes than do one or two measurements made in the vicinity of the tower base. Figure 7 shows a comparison of predicted and observed fluxes, where the observed energy budget has been closed using the modeled G integrated over the EC flux footprint, instead of the observed G; statistical measures of model performance in predicting the flux components H and LE (the components altered in enforcing closure) are listed in Table 5. With this modification the rmsd in predicted H and LE decreased, with the most significant improvement occurring at the ALV2 site where the error in latent heating was reduced to 10% (35 W m−2).
Figure 8a compares flux measurements integrated over the daytime hours (when solar radiation is nonzero) with time-integrated fluxes extrapolated from instantaneous flux predictions assuming a constant daytime evaporative fraction. A statistical description of these comparisons is provided in Table 4. In general, the agreement is good (12% relative error), but slightly poorer than for the instantaneous flux determinations. This may be due in part to diurnal changes in the EF, but also to the fact that biases in eddy covariance measurements tend to amplify when integrated over daily or annual time scales (Moncrieff et al. 1996). As with instantaneous fluxes, agreement in daytime totals improves when the modeled soil heat flux is used to close the observed energy budget (Fig. 8b; Table 5).
b. Spatial and temporal flux variations
While thermal sharpening did not significantly improve statistical comparisons between modeled and measured fluxes, in all cases it did improve model convergence frequency and visual information content in the output flux fields (Fig. 5; nonconvergent pixels appear as isolated or blocks of black cells). Improvement in model convergence is most evident in scenes with strong, narrow discontinuities, such as along the riverbeds in the BESS and ALV2 scenes (Fig. 5a) and the urban street systems in the scenes containing the NORM site (Figs. 5c,d). For the MARE site (Fig. 5c), using bilinearly interpolated thermal inputs resulted in many nonconvergent pixels fringing the high ET river bottoms, where abrupt changes in vegetation-cover fraction are detected at 30-m resolution. Nonconvergence was reduced in employing the standard thermal-sharpening algorithm and was almost eliminated with the convolved sharpening technique.
The enhancement of visual information content attained through thermal sharpening is most pronounced in the vicinity of the NORM supersite, which is located within the city limits of Norman, Oklahoma (Fig. 9). In the 2-week period prior to 10 June 2001, 89 mm of cumulative precipitation was recorded at the NORM weather station, as compared to 6 mm during the 2 weeks before 12 July. The impact of decreased antecedent rainfall on latent heating for 12 July is clear in Fig. 9—only in riparian areas or grounds receiving regular irrigation (golf course, lawns) are high ET rates sustained and in fact enhanced over 10 June rates because of increased atmospheric demand (Fig. 3). With the CS thermal inputs, higher ET from individual fairways in the Westwood Park Golf Course (just south of the airport) can be detected; these fairways would be more heavily irrigated than adjacent areas in the course. Low ET rates are predicted for both days in industrialized areas around the Sooner Fashion Mall (located in the southwest corner, near the cloverleaf interchange) and eastward along the Main Street business corridor.
The residuals to the sharpening polynomial [ΔT̂RAD,low in Eq. (5)] also reveal interesting information regarding spatial anomalies in the scene-average surface temperature–vegetation-cover relationship. Clouds and water bodies will show up as anomalies (both negative—observed surface temperatures are colder than is typical for a given cover fraction). Variations in soil moisture and vegetative stress may also be revealed (positive residuals for soil moisture deficiency and stressed vegetation; negative for wet soils). In the IDAB scene, for example, a backward-L-shaped field with low vegetation cover shows strong negative residuals compared to an adjacent field of similarly low cover (highlighted in Fig. 4). These residuals appear indicative of high soil moisture content, perhaps due to recent irrigation; note that DisALEXI has predicted higher latent heating in this field (Fig. 5b). The extent to which this particular field stands out in the residual map suggests that this may be a valuable tool in identifying soil moisture anomalies over large regions.
A procedure (DisALEXI) for disaggregating large-scale surface flux predictions from the ALEXI model down to much finer spatial scales has been applied to multiscale and multiwavelength remote sensing data collected with the GOES/AVHRR (5 km) and Landsat (30–60 m) satellites. In this modeling system, TIR and vis/NIR wave bands provide information on surface temperature and vegetation-cover status, while kilometer-scale imagery is used to provide atmospheric boundary conditions for meter-scale surface energy balance evaluations. This nested-mapping technique is particularly well suited for operational applications because it does not need local surface meteorological measurements, nor does it require that a wide range in surface conditions be present within the modeling scene, both of which are common limitations among current thermal-mapping approaches.
Disaggregation provides a means for indirectly validating regional-scale flux predictions by breaking them down to scales consistent with the surface source footprint associated with flux-tower observations. When reaggregated using an analytical footprint weighting function, modeled fluxes from DisALEXI agreed with measurements made at seven stations in the OASIS flux network in Oklahoma to within 10% on an instantaneous basis and 15% for daytime totals. In comparison, 5-km instantaneous flux predictions from the regional-scale ALEXI model exhibited a 17% relative error with respect to ground observations, which was strongly dependent on the exact locations of the measurement sites within the model grid cells. Clearly, direct comparison with individual ground-based tower measurements is not an optimal means for assessing the accuracy of regional-scale flux models, especially over complex landscapes.
The utility of an algorithm (DisTrad) for sharpening thermal satellite imagery to the higher resolutions typically associated with vis–NIR bands was also evaluated in the context of disaggregation results. This technique has potential value in that high-resolution thermal data (<1 km) are not currently available with daily coverage and will become even less accessible as the thermal band is being eliminated from future Landsat platforms. Unlike DisALEXI, DisTrad does require a certain amount of variability be present within the scene to develop a decent regression between temperature and vegetation index. Sharpening thermal data inputs to DisALEXI from 60 to 30 m using DisTrad did not significantly improve agreement between modeled and measured fluxes in comparison to inputs processed with a simple bilinear interpolation. This is not surprising, as the typical scale of surface heterogeneity in the vicinity of the flux stations was on the field scale, ∼200–500 m. Greater benefit in terms of observational agreement is to be expected in applying DisTrad to thermal data from MODIS, where the potential for resolution enhancement is 1 km (TIR) to 250 m (vis/NIR); testing of this hypothesis is underway. Nevertheless, the DisTrad sharpening of Landsat thermal data in the current study did result in superior model output in terms of visual information content and model convergence rate, particularly when the residuals to the sharpening function were nominally smoothed through convolution.
There is significant motivation for developing robust and operational nested flux-mapping capabilities. Flux maps at 101–102-m pixel resolution will be useful in precision agriculture, drought and yield forecasting, groundwater modeling, large-eddy simulation, and in detecting functional changes in natural and managed ecosystems in response to climatic and anthropogenic stressors. The scalability in assessment, from local to landscape to regional scales, that nested modeling provides is beneficial in linking cause with effect, which may be occurring on very different spatial scales. Furthermore, methods for upscaling local ground-based observations into a regional context will be critical to the success of ongoing continental- and global-scale flux measurement programs.
Funding for this research was provided primarily by NASA Grant NAG13-99008 and in part by USDA Cooperative Agreement 58-1265-1-043. Suggestions and comments made by three anonymous reviewers greatly improved the clarity and overall presentation of the paper.
Corresponding author address: M.C. Anderson, 1525 Observatory Dr., University of Wisconsin, Madison, WI 53706. Email: firstname.lastname@example.org