Comptes Rendus

Hydrology, environment (Hydrology–hydrogeology)
Spatial distribution of the air temperature in mountainous areas using satellite thermal infra-red data
Comptes Rendus. Géoscience, Volume 343 (2011) no. 1, pp. 32-42.


Understanding the spatial distribution of air temperature in mountainous areas is essential in hydrological modelling. In the Moroccan High-Atlas range, the meteorological stations network is sparse. In order to get additional information, we investigated the thermal infrared data supplied by the Enhanced Thematic Mapper (ETM + ) sensor onboard the Landsat 7 satellite. The brightness temperature derived from ETM+ images is used as a proxy for air temperature to set up a model that describes its spatial distribution. This model accounts for sun location and topographic characteristics derived from the SRTM digital elevation model. It was evaluated on the Rheraya watershed, a 225-km2 region located within the semi-arid High-Atlas mountain range, using two different sources of data. The first data set consists in in-situ air temperature collected by meteorological stations installed during the experiment at various altitudes from 1400 to 3200 m. The second data set is satellite estimates of snow-covered areas (SCA) derived from MODIS images over the whole catchment at 500 m spatial resolution.

Supplementary Materials:
Supplementary material for this article is supplied as a separate file:

Comprendre la distribution spatiale de la température de l’air dans les zones montagneuses est essentielle pour la modélisation hydrologique. Dans la chaîne du Haut-Atlas marocain, la faible densité du réseau de stations météorologiques rend difficile l’obtention de ce paramètre. Dans ce contexte, nous avons utilisé les données infrarouge thermique fournies par le capteur Enhanced Thematic Mapper (ETM +) à bord du satellite Landsat 7. La variation relative de la température de brillance est utilisée pour construire un modèle qui décrit celle de la température de l’air. Le modèle (MSPAT) prend en compte l’angle d’élévation solaire et les caractéristiques topographiques provenant du modèle numérique de terrain SRTM, avec une résolution spatiale de 90 m. MSPAT a été évalué sur le bassin versant de Rheraya, région de 225 km2, située dans le massif semi-aride du Haut-Atlas, en utilisant deux types de données différentes : (1) des chroniques de températures de l’air, mesurées partir de stations météorologiques installées au cours de l’expérience à différentes altitudes de 1400 à 3200 m ; (2) des images de surfaces enneigées (SCA), estimées à partir d’images MODIS acquises à 500 m de résolution spatiale sur l’ensemble du bassin versant.

Compléments :
Des compléments sont fournis pour cet article dans le fichier séparé :

Published online:
DOI: 10.1016/j.crte.2010.11.004
Keywords: Air temperature, Landsat ETM +, Snow-covered areas, Moroccan High-Atlas
Mot clés : Température de l’air, Landsat ETM +, Surface de neige, Haut-Atlas marocain

Abdelghani Boudhar 1; Benoît Duchemin 2; Lahoucine Hanich 1; Gilles Boulet 2; Abdelghani Chehbouni 2

1 Faculté des Sciences et Techniques de Marrakech, avenue A. Khattabi, BP 549, Marrakech, Morocco
2 CESBIO, université de Toulouse, CNRS, CNES, IRD, 18, avenue Edouard-Belin, BPI 280, Toulouse cedex 4, France
     author = {Abdelghani Boudhar and Beno{\^\i}t Duchemin and Lahoucine Hanich and Gilles Boulet and Abdelghani Chehbouni},
     title = {Spatial distribution of the air temperature in mountainous areas using satellite thermal infra-red data},
     journal = {Comptes Rendus. G\'eoscience},
     pages = {32--42},
     publisher = {Elsevier},
     volume = {343},
     number = {1},
     year = {2011},
     doi = {10.1016/j.crte.2010.11.004},
     language = {en},
AU  - Abdelghani Boudhar
AU  - Benoît Duchemin
AU  - Lahoucine Hanich
AU  - Gilles Boulet
AU  - Abdelghani Chehbouni
TI  - Spatial distribution of the air temperature in mountainous areas using satellite thermal infra-red data
JO  - Comptes Rendus. Géoscience
PY  - 2011
SP  - 32
EP  - 42
VL  - 343
IS  - 1
PB  - Elsevier
DO  - 10.1016/j.crte.2010.11.004
LA  - en
ID  - CRGEOS_2011__343_1_32_0
ER  - 
%0 Journal Article
%A Abdelghani Boudhar
%A Benoît Duchemin
%A Lahoucine Hanich
%A Gilles Boulet
%A Abdelghani Chehbouni
%T Spatial distribution of the air temperature in mountainous areas using satellite thermal infra-red data
%J Comptes Rendus. Géoscience
%D 2011
%P 32-42
%V 343
%N 1
%I Elsevier
%R 10.1016/j.crte.2010.11.004
%G en
%F CRGEOS_2011__343_1_32_0
Abdelghani Boudhar; Benoît Duchemin; Lahoucine Hanich; Gilles Boulet; Abdelghani Chehbouni. Spatial distribution of the air temperature in mountainous areas using satellite thermal infra-red data. Comptes Rendus. Géoscience, Volume 343 (2011) no. 1, pp. 32-42. doi : 10.1016/j.crte.2010.11.004. https://comptes-rendus.academie-sciences.fr/geoscience/articles/10.1016/j.crte.2010.11.004/

Version originale du texte intégral

1 Introduction

Air temperature is an important forcing variable for many environmental models, with a large range of applications concerning agronomy, hydrology or ecology. In snow hydrology, air temperature is obviously the main factor to partition liquid from solid precipitation (snowfall). It can be also used as the main driving variable to estimate snowmelt (Bloschl, 1991; Brubaker et al., 1996; Martinec, 1975; Richard and Gratton, 2001). Therefore, it is fundamental to obtain accurate spatially-distributed temperature values for the purpose of hydrological modelling in snowy areas.

Spatial interpolation is widely used to obtain air temperature estimates at any point of interest. Methods can be of deterministic or stochastic nature and differ according to data availability and scale of application. The most common methods are Kriging and Cokriging techniques (Matheron, 1963), Inverse-Distance Weighting (Legates and Willmott, 1990), Splining (Eckstein, 1989; Hutchinson and Gessler, 1994), and Polynomial regression (Myers, 1990). Many studies were carried out to distribute air temperature based on interpolation methods. Over a large study area in Japan, Ishida and Kawashima (1993) report errors between 1 and 3 K for the interpolation of hourly temperature data cokriged with elevation as an external drift factor. Hudson and Wackernagel (1994) also used cokriging with elevation for predicting air temperatures in Scotland, and showed that this technique performed better than simpler interpolation algorithms. Over smaller areas in Sweden, Soederstrom and Magnusson (1995) reported interpolation errors between 0.4 K and 1.6 K. All these studies have pointed out that the results of kriging are highly dependent on the spatial density and distribution of meteorological stations.

In mountainous areas, air temperature is controlled by different factors related to location and topography (Hudson and Wackernagel, 1994). Therefore, the use of interpolation methods cited above may generate substantial errors and biases (Robeson and Willmott, 1993; Willmott et al., 1991). Conventionally, the method of altitudinal lapse rate is used to estimate temperatures at different elevations (Dunn and Colohan, 1999; Singh and Singh, 2001). More sophisticated geostatistical method such as AURELHY (Bénichou and Le Breton, 1987) aims at accounting for the effects of additional topographic factors (slope, aspect). However, such methods require a dense network to be calibrated. This condition is not met in the majority of mountainous areas, especially in developing countries. Meteorological prediction may offer an alternative, but the grid of meteorological models is too coarse (1 to 10 km) for retrieving the spatial pattern of air temperature at a fine resolution over areas with a complex relief.

Satellite thermal infrared data allows retrieving the spatial distribution of the surface skin temperature over extended regions. Several studies, most of them based on coarse resolution data, aimed at using this information to map air temperature at a regional scale. For instance, the daily maximum air temperature field can be retrieved with an accuracy of about 2.5 K using land surface temperatures derived from NOAA-AVHRR (Vogt, 1996; Vogt et al., 1997). Due to the low spatial resolution of these data (1 km), these methods appear difficult to operate in mountainous regions.

In this context, the objective of this article is to investigate thermal infrared satellite data acquired to improve the description of the spatial distribution of air temperature at a fine resolution over mountainous areas. The study area is the Rheraya watershed, High-Atlas mountain range, Morocco. The joint analysis of ETM+ images, meteorological data set and digital elevation model allowed to set up a model (MSPAT) that simulates continuous (daily time step) and spatially distributed (90 m resolution) air temperatures. MSPAT simulations are compared to the altitudinal lapse rate method (LRM). The comparative evaluation of these two models is performed using air temperature recorded from meteorological stations and using snow cover maps derived from MODIS data at a 500m spatial resolution.

2 Datasets

The experiment took place within High-Atlas, a large mountain range in central Morocco (800 km long and 70 km wide). The highest peak is Jbel Toubkal, which culminates at 4174m above sea level. Precipitations in the mountains allow supplying several large irrigated areas in the surrounding plains (Boudhar et al., 2007; Chaponnière et al., 2005; Chehbouni et al., 2008). They are very irregular in time and space, which is typical of semi-arid areas. Snowfalls occur mainly during the winter season (November to April) in the upper parts of watersheds. The contribution of snowmelt to river flows is significant (between 20 and 45% of the total flows according to the characteristics of watershed, Boudhar et al., 2009).

The experimental data set was collected over the Rehraya atlasic sub-catchment (Fig. 1). This head watershed covers a surface area of about 225 km2 and is characterized by a semi-arid and mountainous climate (precipitation of 360 mm/year at the outlet). The elevation ranges from 1084 to 4167 m, and slopes are very steep, with an average grade of 19%.

Fig. 1

The Rheraya watershed and its meteorological station network.

Le bassin versant de Rheraya et le réseau des stations météorologiques.

The experimental data set consist in climatic data, Landat ETM+ and MODIS images, and the SRTM Digital Elevation Model (DEM).

2.1 Climatic data

Air temperatures were recorded by the six meteorological stations located in Fig. 1. Their main topographic characteristics are summarised in Table 1, Supplementary data, online, Appendix A. These characteristics were retrieved from the DEM in the surroundings each meteorological stations. The stations are all located within or in the close vicinity of the Rheraya catchment, with a large range of altitudes (1400 m to 3200 m) and aspects (from 20° northern to 184° southern exposures). The Solar Position Algorithm (Reda and Andreas, 2003) was used to calculate sun position angle.

All these stations were equipped with HMP45C probes (Vaisala, Finland) that measure air temperature and vapour pressure on a half-hourly time step, except at the Caf station where only minimal and maximal temperatures are available. In order to have a maximal sampling, data analysis is performed for daily maximum air temperatures. As illustrated in the Fig. 2, the altitude is a principal factor controlling air temperature variation: an average lapse rate of –0.56 °C per 100 meters was observed during winter 2007/2008.

Fig. 2

Maximal air temperature recorded at the available weather stations in 2007/2008.

Température de l’air maximale, enregistrée dans les stations météorologiques disponibles en 2007/2008.

2.2 Landsat ETM+

The Enhanced Thematic Mapper (ETM + ) sensor is operational onboard the Landsat 7 satellite. ETM+ observes the Earth in seven spectral bands, amongst which one in the thermal infrared region (ETM + 6: 10.40–12.50 μm). ETM+ images are acquired with a regular revisit time of 16 days at a spatial resolution of 30 m in the solar domain and 60 m for the thermal band (Table 2, Supplementary data, online, Appendix A). The scene size is 170 × 183 km when orbiting at an altitude of 705 km. For this study, six images were available over the Rheraya sub-catchment (Table 3, Supplementary data, online, Appendix A).

2.3 MODIS data

The Moderated Resolution Imaging Spectroradiometer (MODIS) is an Earth Observing System (EOS) instrument on board the Terra and Aqua platforms, launched in December 1999 and May 2002, respectively. The sensor observes the Earth from a sun-synchronous position near the polar orbit at an altitude of 705 km. The sensor scans ±55° from nadir in 36 spectral bands in the visible, near- and short-wave and thermal infrared parts of the electromagnetic spectrum. During each scan, 10 along-track detectors per spectral band simultaneously sample the Earth. From its polar orbit, MODIS views the Earth's entire surface ranging from every day at high latitudes to every other day at low latitudes (Justice et al., 1998).

In this article, we used the datasets of the daily surface reflectance product “MOD09GHK”, which is a product computed from the MODIS Level 1B land bands 1-7. Images were ordered from September 2003 to June 2006 through the Earth Observing System (EOS) data gateway (https://wist.echo.nasa.gov). The main characteristics of MODIS images are provided in Table 4, Supplementary data, online, Appendix A. In this study, we used surface reflectances acquired at 500m spatial resolution of 500m in the blue (band 3), the green (band 4) and the short-wave infrared (band 6) part of the electromagnetic spectrum.

2.4 Digital Elevation Model (DEM)

The digital elevation model (DEM) used in this study was derived from the Shuttle Radar Topography Mission (SRTM) data, which consists of a specially modified radar system that flew onboard the Space Shuttle Endeavour during an 11-day mission in February of 2000. It is an international project led by the U.S. National Geospatial-Intelligence Agency (NGA), U.S. National Aeronautics and Space Administration (NASA), the Italian Space Agency (ASI) and the German Aerospace Center (DLR). SRTM obtained elevation data on a near-global scale to generate the most complete high resolution digital topographic database of Earth, including three resolution products, of 1 km and 90 m resolutions for the world, and a 30 m resolution for the US (USGS, 2004). All SRTM data are freely available at: http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp.

In this study, we used the STRM-DEM version 4 acquired at the full 90m spatial resolution over the area covered by landsat scenes (latitudes 30°28′ - 32°24′ North, longitudes 9°14’ and 9°9’ West). Altitude and derivatives (aspect and slope) were calculated through the ENVI DEM extraction module (© 2008 ITT Visual Information Solutions).

3 Methodology

ETM+ brightness temperature was first calculated (Section 3.1), then converted as a proxy of the maximal air temperature images (Section 3.2). The resulting images were analysed together with derivatives of the SRTM digital elevation model for understanding the topographical factors controlling the air temperature variability (Section 3.3). This analysis allowed to set up the MSPAT models that simulates continuous (daily time step) and spatially distributed (90 m resolution) air temperatures (Section 3.4). MSPAT simulations were compared to the altitudinal lapse rate method (LRM, described in Section 3.5). The comparative evaluation of these two models was performed using air temperature recorded from additional meteorological stations and using snow cover maps derived from MODIS data at a 500m spatial resolution. The snow mapping method is described in Section 3.6.

3.1 Brightness temperature calculation

The computation of brightness temperature from ETM + 6 data first requires the conversion of digital numbers to radiances. For this purpose we used the following formula developed by the National Aeronautics and Space Administration (Eq. (1), after Markham and Barker, 1986).

where L(λ) is the spectral radiance received by the sensor (mWcm−2sr−1 μm−1), Qmax is the maximum recorded (255), and Qdn is the digital number of the analysed pixel, Lmin(λ) and Lmax(λ) are the minimum and maximum spectral radiance, detected for Qdn = 0 and Qdn = 255, respectively. For ETM6 with central wavelength of 11.475 mm, it has been set that Lmin(λ) = 0.1238 for Qdn = 0 and Lmax(λ) 1.56 mWcm−2sr−1 μm−1 for Qdn = 255 (Schneider and Mauser, 1996). Thus, the above equation can be simplified into the following form:

Once the spectral radiance L(λ) is computed, the brightness temperature at the satellite level can be directly calculated either by inverting Planck's radiance function for temperature (Sospedra et al., 1998) or by using approximation formula (Eq. (3)), after Goetz et al., 1995; Schott and Volchok, 1985, Wukelic et al., 1989:

where, Tb is the brightness temperature and K1 and K2 are pre-launch calibration constants: K1= 607.76 mWcm−2sr−1 μm−1 and K2 = 1260.56 K, respectively (Schneider and Mauser, 1996).

3.2 Derivation of Maximal Air Temperature from ETM+

Maps of maximal air temperatures were established on the Rheraya watershed by combining ETM+ images and local measurements, with the assumption that air and brightness temperatures (Tb) are linearly related. In principle, there is a tight coupling between the surface and the air temperature under bare soil conditions. Coupling can also be strong when snowpack is below zero degree. Under melting conditions, snow surface is close to zero degree and there is little coupling with the air above. However, Tb distributions according to slope and elevation can be used as additional information to distribute spatially the variations in available melting energy (radiative or sensible heat) in all conditions for use within empirical degree days models1. It is therefore a proxy for these variations rather than a fully deterministic variable. For each image, the brightness temperatures (Tb) of the pixels that include weather stations were extracted and compared with maximal air temperatures (Tam) obtained from the meteorological stations. It appeared that Tb was well correlated with Tam, with a correlation coefficient R2 of 0.82 (Fig. 3). Tam were thus calculated for Tb using a unique regression equation (Tam = 0.7 × Tb–3.2) that was successively applied on each image.

Fig. 3

ETM+6 brightness temperatures versus maximal air temperature measurements.

Température de l’air mesurée en fonction de la température de brillance ETM+6.

3.3 Analysis of Maximal Air Temperature Variation

Fig. 4 represents an example of the co-variation of maximal air temperature and altitude over the Rheraya watershed in 27/02/2003. Not surprisingly, the two variables were found highly correlated, with a decrease by 18 °C of air temperature associated to a 3000 m difference in elevation. This corresponds to an equivalent lapse rate value of 0.60°/100 m. Not shown here, the Tam versus elevation curves extracted for the remaining images display the same pattern, with lapse rates varying between 0.43°/100 m to 0.75°/100 m. These values appear consistent with that of an international standard atmosphere (0.65 °C/100 m according to the International Civil Aviation Organization).

Fig. 4

Air temperature as a function of altitude. The air temperature is calculated from ETM+6 (image of February 27, 2003) and averaged by 5 meters altitudinal zone.

Variation de la température de l’air en fonction de l’altitude. La température de l’air est calculée à partir d’images ETM+6 (image du 27 février 2003) et moyennée par tranche d’altitude de 5 m.

Fig. 5 provides an example of air temperature variations with aspect in February 27, 2003. The aspect values are expressed in degrees, with the values 0° (or 360°), 90°, 180° and 270° corresponding to the facing North, East, South and West, respectively. A difference of 10 °C in air temperature between northern and southern exposures was observed this day. The variations of air temperature with aspect followed a sinusoidal curve centred on about 140°, in coherence with the sun diurnal course (the sun azimuthal angle is 165° at the time of acquisition).

Fig. 5

Aspect effect in air temperature variation. The air temperature is calculated from ETM+6 (image of February 27, 2003) and averaged by 2̊ zone.

Effet de l’exposition des versants sur la variation de la température de l’air. La température de l’air est calculée à partir d’images ETM+6 (image du 27 février 2003) et moyennée par bande de 2.

A deeper investigation shows that the dependence of Tam with aspect is function of altitude, slope and sun elevation. Fig. 6 presents the amplitude of Tam variations with aspects for seven altitudinal zones from 1000 to 4200 m (with steps of 400 m), together with slopes. Variations of Tam with exposures match the altitudinal distribution of slopes, especially during the winter period when the sun elevation is low. For the ETM+ image acquired in May with a maximal sun elevation of 77°, this effect appears very limited.

Fig. 6

Comparison between slope and amplitude of temperature variation for each altitudinal band.

Comparaison entre la pente et l’amplitude de variation de la température de l’air pour chaque bande altitudinale.

3.4 Spatial distribution of Air Temperature (MSPAT model)

The air temperature in any point Tam(i,j) of the Rheraya sub catchment's was calculated from the air temperature recorded at one reference meteorological station (Tam(Iref,Jref), OukaSM and Armed stations, and topographical differences between this point and the reference station. These differences were accounted for using the MSPAT model with the altitude, exposition and slope extracted at the pixel (i,j) and at the location of the reference station (Eq. (4)). The MSPAT was set up to obtain daily images of maximal air temperature at a 90m spatial resolution.


where:G is the air temperature lapse rate calculated from the two reference meteorological stations (Armed and OukaSM) at a daily time step; ΔAlt=Alt(i,j)Alt(ref) is the difference in elevation between the reference station and the pixel where air temperature is calculated.

AS(I,j) and ASref are the aspect in degree of the pixel (i, j) and of the reference station, respectively.

For a given altitude, the model accounts for aspect, slope and sun elevation according to equations (Eqs. (5) and (6)). In details, the air temperature was calculated from the temperature averaged for a given altitudinal band with a correction of the effect of exposure based on a cosine function (Eq. (5)). The maximal air temperature was logically supposed to be the highest when the surface is fully oriented at South (aspect of 180°). The amplitude of the cosine was adjusted as a function of the ratio of the local slope to the sun elevation (A coefficient in Eq. (6) and Fig. 7).

where Tamz is the maximal air temperature for a pixel of altitude z and aspect ASz in degree. Tamz is the mean of maximal air temperature for the altitudinal zone z. A is a coefficient depending of date and altitudinal zone. SE and SL are the sun elevation angle and the slope, both in degree.

Fig. 7

Relationship between the coefficient “A” and the ratio of the sun elevation angle to the average slope of each altitudinal zone for ETM+ date.

Relation entre le coefficient « A » et le rapport angle d’élévation solaire et pente moyenne de chaque tranche altitudinale pour les dates ETM+.

3.5 Lapse rate method (LRM)

In the lowest 10–15 km of the atmosphere, temperature declines with an increase in altitude. This decline, known as temperature lapse rate, is controlled by the balance between heat convection from the surface and radiative cooling. However, lapse rates vary with: (1) altitude, (2) season, (3) latitude, and (4) interactions between topography and weather (e.g. dry versus moist air and inversions).

The lapse rate method (LRM) takes the maximum air temperature observed at a reference station and extrapolates that value over the entire basin via a functional relationship between air temperature and elevation data. This method assumes 1) a linear relationship between air temperature and elevation, and 2) horizontal gradients due to topographic and orographic effects are negligible. In this study, the temperature lapse rates (G) on a daily basis using the stations where the longest dataset is available as:

where Tam1 and Tam2 are the values of maximum air temperature measured at the Armed station and OukaSM stations, respectively. Alt1 and Alt2 are the elevations of these two stations.

Air temperature lapse rates were calculated for each study year and then applied to elevation data for Rheraya basin to determine air temperature for each pixel using air temperature measured at a reference station (Tam(Iref,Jref)) (Eq. (8)).

where: ΔAlt=Alt(i,j)Alt(ref) is the difference in elevation between the reference station (OukaSM and Armed) and the pixel where air temperature is calculated.

3.6 Estimates of snow-covered areas (SCA)

Snow-covered areas were estimated using two methods: (1) Simulation of snowfall and snowmelt, (2) Mapping using MODIS images.

3.6.1 Simulation

The snow cover fraction (SCA) was parameterized as an asymptotically function of the snow water equivalent (SWE) according to equation (Eq. (9)) (Anderson, 1976).

where SWE is the snow water equivalent (in mm), and SCAmax is the maximum allowed snow cover fraction, set to 0.95.

SWE was calculated as the balance between the accumulations processes (snowfall) and the ablation processes (snowmelt). During precipitation events, SWE increased by the daily snowfall if the mean air temperature is below 0 °C, otherwise precipitation is rain. Snowmelt was calculated with a classical temperature index described in Eq. (10).

where Ms is the daily bulk loss of water from the snowpack (mm/day); Td is the mean daily air temperature (°C); T0 is the base temperature (0 °C); and a, the degree-day factor (mm/°C/day). The later is computed as 1.1 multiplied by the ratio of snow and water densities (Martinec, 1975).

Two sets of simulations were performed for years 2003/2004, 2004/2005 and 2005/2006 over the entire Rheraya watershed using a 90 m spatial grid, with the air temperature spatially distributed either with the MSPAT model or with the simple lapse rate method. In both case, the precipitations are distributed using the precipitations recorded at the ARMED station and an average observed elevation gradient of 0.03 mm/100m (Boudhar et al., 2009).

4 Mapping

To select the cloud-free data, we initially eliminated the images showing high reflectances in the blue band (on average higher than 20%) over the lower-lying areas close to the High-Atlas (approximately 5000 snow-free pixels, where the altitude is lower than 1000 meters). Indeed, these high reflectances are due to the presence of clouds on the Atlas foothills, and therefore probably on nearby Atlas summits. In a second phase, we identified clouds on the remaining images by visual examination. After this filtering phase, the number of useful images is respectively 44, 52 and 82 for the seasons 2003/2004, 2004/2005 and 2005/2006 (equivalent of 1 image per week).

The snow-mapping algorithm makes use of the contrast of reflectances between a visible band (MODIS band 4, 0.545–0.565 μm) and shortwave infrared band (MODIS band 6, 1.628–1.652 μm). These two bands were combined to calculate the normalized difference snow index (NDSI, Eq. (11)), which is a proxy of SCA (Hall et al., 2002).

where B4 and B6 are the reflectances in the 4th and the 6th MODIS spectral bands.

In order to reduce errors due to variations of the background reflectance, we used the Modified Normalized Difference Snow Sndex (MNDSI, Eq. (12)). This approach was successfully tested by Chaponniere et al., 2005 using a snow index specifically designed for the SPOT-VEGETATION sensor. This methodology was applied for all cloud-free MODIS images acquired during three years from 2003 to 2005.


Here NDSI0 is the snow index of totally snow-free surfaces (computed, for each pixel, from images acquired in summer), and NDSI100 are the snow indices of pixels totally covered with snow (one single value was fixed as the maximum of NDSI on images acquired during winter).

5 Results and discussion

5.1 Model evaluation at Local scale

Daily maximum air temperatures modelled using MSPAT at the location of meteorological stations were evaluated against all measurements available over the Rheraya watershed, together with estimates of air temperature based on the lapse rate method (referred to as “LRM”). The evaluation was carried out using classical statistical variables (RMSE, BIAS, and R2).

Fig. 8 present the scatter plots between observed and simulated air temperatures during the 2007/2008 season, year with the largest number of available climatic data. Table 5 (Supplementary data, online, see Appendix A) shows the statistical variables obtained by comparing the predictions and the measurements of daily maximum daily air temperature on an annual basis. Biases and RMSE were null for the Armed station which provides with the reference air temperature. Bias was also null for the OukaSM station since it was used for calculating the altitudinal lapse rate. From the statistics calculated on the remaining stations, we noted that the performance of LRM and MSPAT models were similar, with moderate to excellent correlation (R2 from 0.35 to 0.94), low biases (from –0.5° to 1.2°) and error varying from 1 to 3.2°. The MSPAT model performed slightly better than the LRM method for 2 stations (Caf and Neltner), whereas it is the opposite for 2 others (Imsker and Tachdert). These results showed that, at local scale, the MSPAT model did not improve the estimate of air temperature.

Fig. 8

Comparison between observed and simulated air temperature using lapse rate method (LRM) and MSPAT model at different meteorological stations in 2007/2008.

Comparaison entre les températures de l’air observées et simulées, en appliquant la méthode du gradient altitudinal (LRM) et le modèle MSPAT dans les différentes stations météorologiques en 2007/2008.

5.2 Catchment scale

The snow-covered areas derived from MODIS data and simulated with climatic data (LRM and MSPAT models) were compared for the three seasons of study (2003–2004 to 2005–2006). The time series of SCAs averaged on the whole Rheraya watershed were plotted in Fig. 9. The two simulated snow covered area compared favourably to that derived from MODIS images. However, there was a slight underestimation of SCA simulated using the LRM method during the two previous years (2003 to 2005) and a slight overestimation of SCA simulated using the MSPAT model at the beginning of the last year (December 2005). The peaks of SCA were well modelled with the two methods, even though those occurring after generalized snowfall events were clearly underestimated. This smoothing effect was due to the use of an asymptotic function between SCA and snow water equivalent.

Fig. 9

Snow depletion curve estimated with MSPAT model and LRM method in the Rheraya watershed.

Courbe de variation de l’enneigement dans le bassin versant de Rheraya, estimée avec le modèle MSPAT et la méthode LRM.

The Fig. 10 show three examples of SCA maps simulated by the two climatic models and derived from MODIS data during a snowmelt period between 03/12/2003 and 30/12/2003. Simulated SCA appeared much more realistic when the MSPAT model was used to distribute air temperature. This was obvious on the two first dates when looking at the north-eastern and at the western part of the catchment, as well as for the last date for which the spatial extent of snowy areas simulated using the LRM data pattern appeared too limited. This qualitative comparison allowed to show that patterns of SCA were not only matching altitude variations but also exposures. This was consistent with previous studies carried out at the scale of the whole High-Atlas mountains range showing that southern facing slopes have lower SCA than northern ones (Boudhar et al., 2007).

Fig. 10

Snow-covered areas in the Rheraya watershed simulated with the lapse rate method (left, 90 m spatial resolution), MSPAT model (middle, 90 m resolution), and derived from MODIS (right, 500 m resolution).

Surfaces des neiges dans le basin versant de Rheraya, simulées par la méthode du gradient altitudinal (gauche, résolution spatiale de 90 m), le modèle MSPAT (centre, résolution spatiale de 90 m), et dérivées à partir des données MODIS (droite, résolution spatiale de 500m).

The scatter plot between SCA simulated by the two climatic models and derived from MODIS allowed to quantify the previous findings (Fig. 11). This comparison was carried out at a degraded resolution of 2 km2 in order to limit the impact of miss-registration between MODIS and simulated SCA images. SCA modelled with LRM model were found very under-estimated, especially during seasons 2004/2005 and 2005/2006. In contrast, SCA modelled using the MSPAT were less scattered. The statistics results for each season are summarized in Table 6, Supplementary data, online, see Appendix A. We noted a remarkable improvement of the SCA modelled with MSPAT model. This improvement was significant in the two seasons 2004/2005 and 2005/2006 with decreasing in bias by 2.9 and 1.8%, respectively. Theses differences in SCA estimating between LRM and MSPAT models are due principally to the aspect effect not taken account by LRM formulation. In the later, air temperature varies only with altitude and horizontal gradients due to topographic and orographic effects are negligible.

Fig. 11

SCA modelled using lapse rate method and MSPAT method and computed with MODIS data, Comparison pixel by pixel of 2 km2.

Surfaces d’enneigement modélisées en utilisant la méthode du gradient altitudinal et le modèle MSPAT et extraites à partir des données MODIS, comparaison pixel par pixel de 2 km2.

6 Conclusion

This study focused on testing a method of air temperatures spatial variation, principally developed for high-relief mountain environments with scare meteorological data. The thermal infrared band of Enhanced Thematic Mapper (ETM + ) images was used to calculate the brightness temperature. Brightness temperatures were compared with maximal air temperature measured in meteorological stations available, and the correlation equation obtained allowed to compute maps of a proxy of the maximal air temperature. This information was used to characterize its spatial variation according to topographic factors (elevation, aspect and slope) and solar elevation angle. From this, we set up a model to spatialise air temperature.

The results of air temperature simulation with MSPAT model was compared with lapse rate model (LRM). The simulated air temperatures with the two models were similar at the location of available meteorological stations. Over all the study area, the models were evaluated using snow cover area (SCA) estimates. The last are then compared with SCA product mapped from MODIS data. Results showed that MSPAT model had a considerable improvement on the modelling SCA. It is thus believed that the pseudo air temperature approximated using ETM+ brightness temperature may be useful to monitor the dynamics of snowpacks. The approach developed here may be beneficial for hydrological applications, especially snowmelt modelling at watershed scale with scare data.


This study was supported by the research projects SUDMED (IRD-UCAM), PAI (Programme d’Action Intégrée du Comité Mixte Interuniversitaire Franco–Marocain, Jeune Equipe IRD, CREMAS), Volubilis (MA/06/148) and the PLEADeS Project (http://www.pleiades.es/) funded by the European Commission (Sixth PCRD).

1 Boudhar, A., Boulet, G., Hanich, L., Sicart, J.E. and Chehbouni, A., submitted. Measurement and modeling of the energy fluxes and melt rate of a seasonal snow cover at the Oukaimden site in the Moroccan High Atlas. Hydrology and Earth System Sciences.


[Anderson, 1976] Anderson, E.-A., 1976. A point energy and mass balance model of a snow cover. Silver Spring, MD US. National Oceanic and Atmospheric Administration (NOAA), Technical Report NWS 19.

[Bénichou and Le Breton, 1987] P. Bénichou; O. Le Breton Prise en compte de la topographie pour la cartographie des champs pluviométriques statistiques, La Météorologie, 1987 (14 p)

[Bloschl, 1991] G. Bloschl The influence of uncertainty in the air temperature and albedo on snowmelt, Nordic Hydrol., Volume 22 (1991), pp. 95-108

[Boudhar et al., 2007] A. Boudhar; B. Duchemin; L. Hanich; A. Chaponnière; P. Maisongrande; G. Boulet; J. Stitou; A. Chehbouni Snow covers dynamics analysis in the Moroccan High-Atlas using SPOT-VEGETATION data, Sécheresse, Volume 18 (2007) no. 4, pp. 278-288

[Boudhar et al., 2009] A. Boudhar; L. Hanich; G. Boulet; B. Duchemin; B. Berjamy; A. Chehbouni Impact of the snow cover estimation method on the Snowmelt Runoff Model performance in the Moroccan High-Atlas Mountains, Hydrol. Sci. J./J. Sci. Hydrol., Volume 54 (2009) no. 6, p. 2009

[Brubaker et al., 1996] K. Brubaker; A. Rango; W. Kustas Incorporating radiation inputs into the snowmelt runoff model, Hydrol. Process, Volume 10 (1996), pp. 1329-1343

[Chaponnière et al., 2005] A. Chaponnière; P. Maisongrande; B. Duchemin; L. Hanich; G. Boulet; S. Halouat; R. Escadafal A combined high and low spatial resolution approach for mapping snow covered area in the Atlas Mountain, Int. J. Remote Sensing, Volume 26 (2005), pp. 2755-2777

[Chehbouni et al., 2008] A. Chehbouni; R. Escadafal; G. Boulet; B. Duchemin; V. Simonneaux; G. Dedieu; B. Mougenot; S. Khabba; H. Kharrou; O. Merlin; A. Chaponnière; J. Ezzahar; S. Erraki; J. Hoedjes; R. Hadria; H. Abourida; A. Cheggour; F. Raibi; A. Boudhar; L. Hanich; N. Guemouria; Ah. chehbouni; A. Olioso; F. Jacob; J. Sobrino An integrated modelling and remote sensing approach for hydrological study in semi-arid regions: the SUDMED Program, Int. J. Remote Sensing, Volume 29 (2008), pp. 5161-5181

[Dunn and Colohan, 1999] S.M. Dunn; R.J.E. Colohan Developing the snow component of a distributed hydrological model: a step-wise approach based on multiobjective analysis, J. Hydrol., Volume 223 (1999), pp. 1-16

[Eckstein, 1989] B.A. Eckstein Evaluation of spline and weighted average interpolation algorithms, Comput. Geosci., Volume 15 (1989) no. 1, pp. 79-94

[Goetz et al., 1995] S.J. Goetz; R.N. Halthore; F.G. Hall; B.L. Markham Surface temperature retrieval in a temperate grassland with multiresolution sensors, J. Geophys. Res., Volume 100 (25) (1995) (397-25 410)

[Hall et al., 2002] D.K. Hall; G.A. Riggs; V.V. Salomonson; N.E. DiGirolamo; K.J. Bayr MODIS snow-cover products, Remote Sensing Environ., Volume 83 (2002), pp. 181-194

[Hudson and Wackernagel, 1994] G. Hudson; H. Wackernagel Mapping temperature using kriging with external drift: theory and an example from Scotland, Int. J. Climatol., Volume 14 (1994), pp. 77-91

[Hutchinson and Gessler, 1994] M.F. Hutchinson; P.E. Gessler Splines – more than just a smooth interpolator, Geoderma, Volume 62 (1994), pp. 45-67

[Ishida and Kawashima, 1993] T. Ishida; S. Kawashima Use of cokriging to estimate surface air temperature from elevation, Theor. Appl. Climatol., Volume 47 (1993), pp. 147-157

[Justice et al., 1998] C. Justice; E. Vermote; J.R.G. Townshend; R. Defries; D.P. Roy; D.K. Hall; V.V. Salomonson; J. Privette; G. Riggs; A. Strahler; W. Lucht; R. Myneni; Y. Knjazihhin; S. Running; R. Nemani; Z. Wan; A. Huete; W. van Leeuwen; R. Wolfe; L. Giglio; J.-P. Muller; P. Lewis; M. Barnsley The Moderate Resolution Imaging Spectroradiometer (MODIS): land remote sensing for global change research, IEEE Trans. Geosci. Remote Sensing, Volume 36 (1998), pp. 1228-1249

[Legates and Willmott, 1990] D.R. Legates; C.J. Willmott Mean seasonal and spatial variability in global surface air temperature, Theor. Appl. Climatol., Volume 41 (1990), pp. 11-21

[Markham and Barker, 1986] Markham, B. L., and Barker, J. L., 1986. Landsat-MSS and TM post calibration dynamic ranges, atmospheric reflectance and at-satellite temperature. EOSAT Landsat Technical Notes 1, (Lanham, Maryland: Earth Observation Satellite Company), pp. 3–8.

[Martinec, 1975] J. Martinec Snowmelt-runoff model for stream flow forecasts, Nordic Hydrol., Volume 6 (1975), pp. 145-154

[Matheron, 1963] G. Matheron Principles of geostatistics, Econ. Geol., Volume 58 (1963), pp. 1246-1266

[Myers, 1990] R.H. Myers Classical and Modern Regression with Applications, PWS-Kent Publishing, Boston, MA, 1990 (488 p)

[Reda and Andreas, 2003] Reda, I., Andreas, A., 2003. Solar position algorithm for solar radiation application. National Renewable Energy Laboratory (NREL) Technical report NREL/TP-560-34302.

[Richard and Gratton, 2001] C. Richard; D.J. Gratton The importance of the air temperature variable for the snowmelt runoff modelling using the SRM, Hydrol. Process., Volume 15 (2001), pp. 3357-3370

[Robeson and Willmott, 1993] Robeson, S. M. and Willmott, C. J., 1993. Spherical spatial interpolation and terrestrial air temperature variability, Proceedings Second International Conference on Integrating GIS and Environmental Modeling, Breckenridge, CO, pp. 111–115.

[Schneider and Mauser, 1996] K. Schneider; W. Mauser Processing and accuracy of Landsat ThematicMapper data for lake surface temperature measurement, Int. J. Remote Sensing, Volume 17 (1996), pp. 2027-2041

[Schott and Volchok, 1985] J.R. Schott; W.J. Volchok Thematic Mapper thermal infrared calibration, Photogramm. Eng. Remote Sensing, Volume 51 (1985), pp. 1351-1357

[Singh and Singh, 2001] P. Singh; V.P. Singh Snow and glacier hydrology, Klüwer, Dordrecht, 2001

[Sospedra et al., 1998] F. Sospedra; V. Caselles; E. Valor Effective wave number for thermal infrared bands – application to Landsat-TM, Int. J. Remote Sensing, Volume 19 (1998), pp. 2105-2117

[USGS, 2004] USGS, 2004. Shuttle Radar Topography Mission, 3 Arc Second scene SRTM_u03, Unfilled Unfinished 2.0, Global Land Cover Facility, University of Maryland, College Park, Maryland, February 2000.

[Vogt, 1996] Vogt, J.V., 1996. Land surface temperature retrieval from NOAA AVHRR data. In: D'Souza, G., Belward, A. & Malingreau, J.P. (Eds.). Brussels & Luxembourg, pp. 125–151.

[Vogt et al., 1997] J.V. Vogt; A.A. Viau; F. Paquet Mapping regional air temperature fields using satellite derived surface skin temperatures, Int. J. Climatol., Volume 17 (1997), pp. 1559-1579

[Willmott et al., 1991] C.J. Willmott; S.M. Robeson; J.J. Feddema Influence of spatially variable instrument- networks on climatic averages, Geophys. Res. Lett., Volume 18 (1991), pp. 2249-2251

[Wukelic et al., 1989] G.E. Wukelic; D.E. Gibbons; L.M. Martucci; H.P. Foote Radiometric calibration of Landsat Thematic Mapper thermal band, Remote Sensing Environ., Volume 28 (1989), pp. 339-347

Comments - Policy