A Wetness Index Using Terrain-Corrected Surface Temperature and Normalized Difference Vegetation Index Derived from Standard MODIS Products: An Evaluation of Its Use in a Humid Forest-Dominated Region of Eastern Canada.

In this paper we develop a method to estimate land-surface water content in a mostly forest-dominated (humid) and topographically-varied region of eastern Canada. The approach is centered on a temperature-vegetation wetness index (TVWI) that uses standard 8-day MODIS-based image composites of land surface temperature (TS) and surface reflectance as primary input. In an attempt to improve estimates of TVWI in high elevation areas, terrain-induced variations in TS are removed by applying grid, digital elevation model-based calculations of vertical atmospheric pressure to calculations of surface potential temperature (θS). Here, θS corrects TS to the temperature value to what it would be at mean sea level (i.e., ∼101.3 kPa) in a neutral atmosphere. The vegetation component of the TVWI uses 8-day composites of surface reflectance in the calculation of normalized difference vegetation index (NDVI) values. TVWI and corresponding wet and dry edges are based on an interpretation of scatterplots generated by plotting θS as a function of NDVI. A comparison of spatially-averaged field measurements of volumetric soil water content (VSWC) and TVWI for the 2003-2005 period revealed that variation with time to both was similar in magnitudes. Growing season, point mean measurements of VSWC and TVWI were 31.0% and 28.8% for 2003, 28.6% and 29.4% for 2004, and 40.0% and 38.4% for 2005, respectively. An evaluation of the long-term spatial distribution of land-surface wetness generated with the new θS-NDVI function and a process-based model of soil water content showed a strong relationship (i.e., r2 = 95.7%).


Introduction
Soil water content (SWC) is a measure of the total amount of water, including water vapor, contained within a column of soil above the ground water table. Water diffuses into the soil through the soil pores and is transferred through the soil at a rate defined by the soil's permeability. Available soil water is used by plant growth processes and nutrient acquisition food production (along with available light energy) and in many instances is released directly to the atmosphere through the plant-mediated processes of transpiration and evaporation of free water in the soil. SWC is a critical site variable that defines among many factors (a) the duration and intensity of drought, (b) agricultural and tree crop production and food security, (c) level of soil erosion and surface runoff impacts, and (d) rate of evapotranspiration (ET). In particular to forest-dominated regions, SWC plays a vital role in defining (a) the condition and health of forests, (b) the distribution and rate of forest growth and regeneration, (c) forest fire hazard rating, (d) level of carbon dioxide assimilation by forests, and (e) scheduling of forest harvesting and other forest management activities, such as aerial seeding and tree planting.
The most commonly used and standard protocols for estimating SWC are the "gravimetric" and "volumetric" methods. Both methods estimate SWC accurately, but only provide point measurements due to the amount of work required to process a single field sample. Automated point-sampling possible with time domain reflectometry and other electronic devices can improve the spatial coverage, but at an enormous cost. In general, field measurements of SWC are insufficient to provide the near continuous spatial coverage needed for land-management applications. As an alternative to in-field measurement of SWC, remote sensing (RS) platforms, like the TERRA satellite, can be used to generate a vast amount of information about the earth surface at resolutions small enough and frequent enough (everyday with the MODIS sensor) to be appropriate for land-management planning.
For the last two decades, optical and thermal RS sensors have been employed for estimating water content (WC) in soils and vegetation. An overview by Moran et al. [1] reveals that a strong relationship exists between WC and index values derived from surface temperature-vegetation index (T S -VI) methods that use optical and thermal RS data as input. Although showing promise, the penetrating capabilities of optical and thermal RS sensors to sample the earth's surface diminish with the presence of clouds, vegetation cover, and other non-transmitting surfaces. At present, many challenges remain with employing RS data in understanding earth-system processes.
The relation between T S and VI has been previously described in scatterplots of T S vs. VI, typically generating either "triangular" [2][3] or "trapezoidal" forms [4][5]. The edges of these forms provide the basis for modelling WC in bare soils and in vegetation (see Section 3.2 for more detail). A negative slope of the T S -VI plot is generally observed and interpreted using micro-meteorological principles. For instance, low T S 's generally occur over dense vegetation with high SWC, and relatively higher T S 's, over sparse vegetation with low SWC. In the literature, T S -VI plots have been given many different names, such as the surface moisture status [6], triangular method of indirect estimation of soil moisture [7], surface temperature-vegetation index [5], moisture index (MI; [8]), temperature-vegetation contextual (TVX; [9]), and temperature-vegetation dryness index (TVDI; [3]). The T S -VI approach has been implemented using various RS data over a diversity of ecosystems worldwide ( Table 1). The approach has also been applied in various applications, such as in quantifying (i) the evaporative fraction [10][11]), (ii) ET [2,12]), (iii) crop water deficit [4], (iv) forest fire hazard [13], and (v) drought duration and intensity [14].
The widely used T S -NDVI relationship has limited application in variable topography [17]. For instance, remotely-sensed T S can only be estimated under clear-sky conditions and T S in variable terrain is often lower in upland areas compared to low-lying areas within the same climatic and landuse zonation. Clearly, using a response function based on strictly T S could incorrectly suggest that land-surface conditions in upland areas are wetter than in low-lying areas. In an effort to improve WC calculations in high elevation areas, corrections are applied to T S by applying grid, digital elevation model (DEM)-based calculations of vertical atmospheric pressure to calculations of surface potential temperature [θ S ; 18]. Here, θ S adjusts T S to a value representing the temperature that would be at mean sea level (i.e., ~101.3 kPa) in a neutrally-stratified atmosphere [18].
The objectives of this paper are to: (i) develop a temperature-vegetation wetness index (TVWI) based on an application of an θ S -NDVI response function constructed from MODIS-based standard products of T S and surface reflectances; and (ii) evaluate the potential of TVWI to estimate both the spatial and temporal variations in SWC in a humid forest-dominated and topographically-variable landscape in eastern Canada.

Ref.
T S -VI approaches used in the past * [3] • Assumed variations in T S were caused by variations in SWC; TVDI was derived from a T S -NDVI * relation using NOAA-AVHRR images taken over Senegal in West Africa.
• Generally a good agreement existed between estimates of SWC from the T S -NDVI function and estimates from hydrological models using ground-based climate data as input to the models (r 2 =70%). [5] • Based on a trapezoidal form of the T S -F R * function using NOAA-AVHRR images taken over Africa.
• The method illustrated the conceptual borders of dry and wet soils and areas of low-to-high ET with respect to T S and vegetation coverage. [6] • Examined the effect of biome type on the slope of the T S -NDVI function at regional and continental scales using NOAA-AVHRR images taken over western Montana, USA. • The approach provided a strong relationship between the slope of the T S -NDVI function and a crop moisture index (yielding an r 2 of 83%), pointing to the fact that the T S -NDVI function is sensitive to surface soil water conditions. [7,15] • Represented SWC using a "triangular" approach, based on a function of T S and NDVI using airborne images taken over the states of Kansas and Arizona, USA and processed with the SVAT * model. • Method provided a standard error of 16% between T S -NDVI wetness-value predictions and ground measurements of SWC. [8] • Used T S and a ratio of near infrared and the blue bands for modelling MI using LANDSAT images taken over northeast Brazil.
• Method modelled the seasonal variations of SWC due to changes in vegetation growth; correlations between MI and ground-based measurements of SWC were not reported. [9] • Established relationships between the slopes of T S -NDVI (i.e., TVX approach) and soil moisture produced with the Simplified Simple Biosphere Model using NOAA-AVHRR images taken over the HAPEX-Mobilhy site in southern France.
• TVX metrics provided estimates of daily soil moisture variation up to 2-cm depth in the soil and seasonal variations up to 10 cm.
• Method provided a general linear relation with the residual standard error of < 4% over a range of 6-27% in volumetric SWC. [16] • Modelled SWC from a T S -F R relation using NOAA-AVHRR and LANDSAT ETM+ images taken over northeastern Spain.
• Method provided a reasonable agreement between NOAA and LANDSAT-based estimates of SWC (r 2 =66%).

Study Area and Data Used
Canada is divided into fifteen terrestrial ecozones or generalized zonal groupings of similar soil formation, climate, and landuse cover [19]. In eastern Canada, the Atlantic Maritime Ecozone is dominated by forests occupying approximately three-fourths of the total ecozone's land-surface area. The combined surface area of the Provinces of New Brunswick (NB), Nova Scotia (NS), and Prince Edward Island (PEI) occupy about 80% of the total landbase of the Atlantic Maritime Ecozone. The forest in the region is generally characterized by a temperate mix wood forest of deciduous species, such as maple (Aceraceae spp.), beech (Fagus grandifolia Ehrh.) and birch (Betulaceae spp.), and coniferous species, such as spruce (Picea spp.) and balsam fir [Abies balsamea (L.) Mill (1) the westerly displacement of continental air masses, (2) the region's proximity to the Atlantic Ocean, and in the case of NB's southern shore and NS' northern shore, by the cold water of the Bay of Fundy, and (3) changes in elevation. The area experiences a cool-moist climate with a mean annual temperature and annual precipitation range of 3.5-6.5°C and 900-1500 mm, respectively. The Provinces of NB and NS have elevational ranges from 0-834 m for NB, and 0-540 m above mean sea level (a.s.l.) for NS. PEI is relatively flat, having an average elevation of about 30 m a.s.l. (Figure 1).
In this study, we employed 8-day composites of MODIS-based standard products of (i) daytime T S at 1-km resolution (MOD11A2 products) and (ii) red and near infrared (NIR) surface reflectances at 250-m resolution (SR; MOD09Q1 products) obtained from NASA to generate products of NDVI at 250-m resolution. Generally, in humid, forest-dominated regions of the world day-to-day variation in SWC seldom control the long-term processes involved with defining forest productivity. As a result, WC values in this paper were calculated over 8-day intervals roughly to coincide with the 8-day standard image composites of T S and SR. To calculate vertical atmospheric pressure variation and θ S [18], we employed a 250-m resolution DEM of the study area generated from 3-arc second resolution height point-data acquired from the NASA Shuttle Radar Topography Mission (SRTM) archive. For comparing MODISderived wetness values, field measurements of volumetric soil water content (VSWC) were acquired from ten sites in southwest NB, including two Fluxnet-Canada sites [20], namely the Nashwaak Lake (NWL) and Charlie Lake  Figure 1). Measurements of VSWC were obtained at all sites with water content reflectometers (i.e., CS615 or CS616 probes, Campbell Scientific Corp.) placed at 10-cm depths. To coincide with the 8-day averaging scheme employed earlier, all field measurements from the individual sites were temporally-averaged to 8-day means. Table 2 provides the time periods over which individual 8-day mean calculations of TVWI and VSWC apply. As a second comparison, we generated spatially-explicit long-term average calculations of SWC (given in "% of saturation") at 220-m resolution for the Province of NB with the well-established, process-based WET model [21][22]. Model output was calculated based on inputs of long-term average values of precipitation (mm d -1 ), soil infiltration capacity (mm d -1 ), solar energy input and exchange (and, therefore, DEM; MJ d -1 ) [23], flow accumulation (in m 2 or number of cells), and surface run-off rates (mm d -1 ) [21][22].   Figure 2 shows a schematic diagram of the TVWI method. The method consists of (1) data preprocessing in the generation of the 8-day composites of NDVI and in the correction of T S -and NDVIimage products by gap filling of cloud-contaminated pixels, (2) calculating θ S and TVWI for 58 images, and (3) performing a comparison of TVWI with field measurements of VSWC and model calculations of SWC ("% of saturation").

Pre-processing of the MODIS Images
Eight-day composites of NDVI were determined using: where ρ is the SR values for the red and NIR bands. Cloud-contamination of pixels was identified by quality control flags contained in MOD09Q1. The 8-day datasets of T S had clear distinction between clear-sky and cloud-contaminated pixels. To correct T S and NDVI values for cloud-contaminated pixels, we employed an image-correction algorithm first proposed by Hassan et al. for the correction of T S [24]. Mathematically, this correction is expressed as where ) i ( X is the mean value of either T S or NDVI for each of their respective 8-day composites for a specific year (where i=1…n, and n=the total number of 8-day composite images involved per year), X(i) represents the timeseries of image-based T S and NDVI values for a specific year for cloud-contaminated pixels taking into account only T S and NDVI values from cloud-free composites, m is the total number of cloud-free pixels from individual composite images available for a specific year for a cloud-contaminated pixel, A is the average temporal deviation of T S and NDVI from their respective ) i ( X 's for a specific cloud-contaminated pixel, and B n is the estimated 8-day means of either T S or NDVI for individual cloudcontaminated pixels.

Calculating θ S and TVWI
There were two steps involved in the calculation of θ S . The first step involved the calculation of atmospheric pressure (p) at each of the image pixels by referring to their respective point elevations in the acquired DEM of the study area and 26 . 5 [25], where p (in kPa) is the atmospheric pressure and z (in m) is the elevation a.s.l. This equation is based on a simplified form of the ideal gas law for a neutrally-stratified atmosphere and a temperature of 20 o C at a standard atmosphere (i.e., 101.3 kPa). Surface potential temperature, θ S (in K), was then calculated with where T S is the image-pixel surface temperature (in K), p o is the average pressure at mean sea level (set at 101.3 kPa), R is the gas constant (i.e. 287 J kg -1 K -1 ), and c p is the specific heat capacity of air (~1004 J kg -1 K -1 ). In normal practice, meteorologists set p o to 100.0 kPa instead of 101.3 kPa that we used. The meteorological principles behind Eq. (4), however, remain the same. As a result of corrections to T S (yielding θ S ), our TVWI is based on a response function having θ S and NDVI as the main axes. A trapezoidal shape (see Figure 3) was generally observed when θ S and NDVI data were plotted (see Figure 5, Results & Discussion section). In Figure 3, the dry edge (θ dry ; the line where θ S is highest in relation to NDVI) represents the case where water is not available for ET, and, as a result, TVWI possesses the lowest value (~0.0). The dry edge (θ dry ) is determined as a linear fit of the highest θ S in relation to NDVI. In contrast, the wet edge (θ wet ; the line where θ S is lowest in relation to NDVI; Figure 3) represents the case where water is freely available for ET and TVWI is highest (~1.0). Pixel-calculations of TVWI are based on a modified form of the wetness index proposed in [3], i.e., wet dry where TVWI is a dimensionless quantity and θ s , θ dry and θ wet are all in K.

Figure 3.
Trapezoidal form of the TVWI based on a relationship between θ S and NDVI (modified after [5]). Dry and wet edges of the trapezoid are indicated.

TVWI and SWC comparisons
There were 58 TVWI maps generated for the years 2003-2005. For comparison of the calculated TVWI at the field level, we determined the point-location estimates of TVWI by (i) taking into account all TVWI-values within a 1 km ×1 km area (representing a 4×4 pixel sampling window) centered on the measurement sites, and (ii) averaging their respective values (16 in all). We analyzed the trend-lines for both the TVWI and VSWC as a function of time period (Table 2) to establish whether temporal trends between the two quantities were similar.
As a second comparison, a map of long-term averages of TVWI was created by averaging the 58 TVWI images. These values were then compared with SWC values produced with the WET model. The WET model gives SWC in "% of saturation" in the range of 0-1, with 0 to 1 representing excessively dry soil conditions to progressively wetter conditions. For each of the "% of saturation" values, we computed the mean TVWI-values from the long-term averages. Least squares regression and the coefficient of determination (r 2 ) were used to test the degree of agreement between the model-derived values of SWC and TVWI. No evaporation θ dry1 = a 1 θ dry2 = a 2 +b 1 *NDVI slight over-prediction with NDVI-values < 0.75 and an under-prediction with NDVI-values > 0.75. In a similar study of T S , an appreciably better correlation was established; yielding r 2 =88.3% and regression slope of 0.86 [24]. A possible reason for this difference may be related to the tendency of NDVI to saturate near the middle of the growing season (Figure 4b), while mean temperatures continue to vary [24]. At the point of saturation, NDVI becomes decoupled from temperature.  Figure 5 shows some example scatterplots (and related dry and wet edges) obtained by plotting θ S against NDVI. Among the cases considered, all scatterplots matched Figure 3, although slight differences in dry and wet edges were observed. Despite θ wet varied among periods (with a range of 268-282 K), we used a constant mean value of 275 K for θ wet to calculate variations in TVWI. The rationale behind this strategy was to ensure that: (i) the same baseline was used to find the maximum TVWI (i.e., 1.0) across the entire series of images, and (ii) the value of θ wet ≥ 275 was sufficiently high to support ET.  of related dry and wet edges (following Figure 3). Figure 6 shows the spatiotemporal variability of VSWC and derived TVWI specifically for the field sites (i.e., at NWL, CL, and CFB Gagetown in Figure 1). For each period of 2003-2005 (Table 2), we generated a pair of box-plots to demonstrate the spatial variability among sites. In general, a larger spatial variability was observed in VSWC compared to TVWI. This was expected given VSWC was acquired at point locations, while TVWI values were obtained by averaging individual pixel values within an 1 km × 1 km area centered over the point of measurement. The mean values represented by the red and blue closed symbols in Figure 6, however, were consistently close (within ~ ±25% or better).  affected period in 2003, we found T S to suddenly drop, which caused TVWI to rise above normal conditions and become somewhat unrepresentative. In 2005, the 8-day T S composite for that period over the field sites ( Figure 1) were generated using one to three days of T S data as a result of excessive cloudiness, which may have provided an unrepresentative calculation of T S and, therefore, TVWI. Apart from these two apparent anomalies, our analysis revealed that the temporal trend lines for VSWC and TVWI were fairly close to each other (Figure 7). The y-intercepts of VSWC and TVWI trend lines indicated a fairly consistent relationship between mean quantities, i.e., 31.0% and 29.0% for 2003, 28.   Table 2, for actual dates). Circled symbols in panels (a) and (c) represent anomalous data removed from trend analysis (see main text for details).  From this analysis, TVWI (on the basis of relations between θ S and NDVI) could be perceived as an indirect measure of SWC because of the following reasons:

Results and Discussion
The TVWI values were extremely close to VSWC values in terms of their magnitudes, even in heavily-forested areas (see Figure 7).
As TVWI is based on the measure of surface properties, TVWI could be viewed as an indirect indicator of canopy (overstory foliage) WC. Canopy foliage WC in normal conditions is in equilibrium with SWC due to the relation between the water balance and total leaf area index [26]. This equilibrium could proceed uninterrupted if plant growth is unaffected by disturbance agents such as disease, insect infestation, physiological drought caused by soil-water logging or disruption of root to shoot ratios by air pollutant deposition [27][28] or winter conditions such as prolonged winter thaws [29]. Moreover, plant water potential (a measure of in-canopy WC) and soil water potential (an indirect measure of SWC) were shown to be strongly correlated [30][31].
In humid, water surplus conditions it is expected that ET would proceed at sufficiently elevated rates to cause measurable cooling of the immediate environment; on average, ~ 2/3 of available net radiation goes to ET in water surplus areas of the world [32]. Figure 9a provides a map of long-term averages of TVWI produced by averaging the 58 TVWI images. TVWI values fell mostly in the range of 20-60% with an average of 40%. In summary: The areas along the coastlines and wider river channels showed higher TVWI, > 50%, due to their proximity to water and low elevation relative to exposed water. (ii) Some high relief areas, such as in northwestern NB and in the eastern part of NS, had comparatively lower TVWI values than other regions in the Maritime Provinces (20-32%). A similar pattern was observed with "% of saturation" values generated with the WET model for the Province of NB, despite some minor differences in low-lying areas of the Province ( Figure  9c). (iii) A comparison of mean values of TVWI with values of "% of saturation" for the same pixels and region provided reasonable agreement (r 2 =95.7%; Figure 10a); indicating that for the most part (for > 99.5% of the values) the wetness distributions in Figures 9a and 9c for the Province of NB, despite having different wetness units (Figures 9b and 10b), were mostly the same. Values of "% of saturation" < 0.45 were not included in the comparison as (i) the mean TVWI values did not show any clear pattern with modelled wetness values (Figure 10a), and (ii) the amount of data involved was fairly small, comprising < 0.5% of total available data points (Figure 10b). Reasons for divergence in the low "% of saturation" range ( Figure 10a) remain unknown to us.  The extended left-hand tail of the frequency distribution in (b) and boxed data points in (a) accounts for < 0.5% of all available data points.

Concluding Remarks
In this paper, we demonstrated a practical way to estimate the spatiotemporal variability in surface wetness conditions in a forest-dominated, topographically-diverse region of eastern Canada using MODIS data. Here, the wetness index (TVWI) is based on a function of θ S -NDVI obtained by processing standard MODIS-based products of 8-day composites of surface temperature and surface reflectance. The use of θ S in the calculation of TVWI provides a novel approach to dealing with terrain-induced variation in temperature and, therefore, TVWI and surface-wetness mapping from satellite optical and thermal imaging of forested landscapes. In order to calculate θ S , we apply a correction to surface temperatures by employing basic principles of meteorology. The TVWI developed here was shown to capture a significant portion of (i) the spatiotemporal variability in field measurements of SWC at ten, mostly forested sites in southwest NB, and (ii) the spatial variability in modelled values of "% of saturation" generated for the Province of NB. The TVWI can be viewed as an indirect measure of SWC as surface vegetation water content in normal, humid conditions has been previously shown by others to be mostly in equilibrium with SWC. point data free of charge. We would also like to thank Mr. Don Coleman, Meteorological Service of Canada, for providing measurements of SWC at field sites managed by CFB Gagetown personnel. Partial financial support was provided to QKH through the "Dorothy and Kenneth Langmaid Forest Soils and Tree Growth Scholarship" administered by the Fundy Community Foundation of NB, Canada.