Interactive comment on “ Combining remote sensing and GIS climate modelling to estimate daily forest evapotranspiration in a Mediterranean mountain area ”

Authors reply to reviewer general comments 1: In the case of the B parameter, we have chosen three ways to compute it depending on its complexity. We know that B parameter computed using a net radiation sensor needs more intensive local calibration. But, on the other hand, B computed using NDVI do not require local calibration, only realistic values of NDVI. When we compared the different B parameter methodologies, results, in Landsat case, were so close, so we can use B-NDVI in an operative way without local calibration. However, in this article we also wanted to point out not only the advantages of the methodology but also the limitations of these techniques in order to improve them in a future.


Introduction
Evaporation and transpiration are the two main processes involved in water transfer from vegetated areas to the atmosphere.Evapotranspiration (ET) from the Earth's vegetation constitutes 88 % of the total terrestrial evapotranspiration, and returns more than 50 % of terrestrial precipitation to the atmosphere (Oki and Kanae, 2006); therefore it plays a key role in both the hydrological cycle and the energy balance of the land surface.Climate warming may accelerate the hydrological cycle as a result of enhanced evaporative demand in some regions where water is not limiting (Jung et al., 2010).However, the combination of warmer temperatures with constant or reduced precipitation in other regions may lead to a large decrease in water availability for natural and agricultural systems as well as for human needs (Jackson et al., 2001), especially in arid or semiarid areas (Jung et al., 2010) such as the Mediterranean basin (Bates et al., 2008).
Evapotranspiration has been measured extensively at local scales using micrometeorological (such as eddy-covariance or the Bowen ratio) or sap flow techniques.Since the last decade, there have been several global initiatives to monitor evapotranspiration in different vegetation types, such as FLUXNET (ORNL DAAC, 2010).Therefore, magnitudes and controls (climate, water availability, physiological regulation, etc.) on evapotranspiration are widely known for different types of vegetation, albeit at small spatial scales.How can we improve, then, our knowledge of evapotranspiration?In terms of spatial variability (and its driving factors) the next challenge is larger scales.
ET can be modelled at global scales using GIS climatebased methodologies such as GIS-based Environmental Policy Integrated Climate -GEPIC - (Liu et al., 2007(Liu et al., , 2009)), Lund-Potsdam-Jena managed Land -LPJmL - (Rost et al., 2008) or Global Crop Water Model -GCWM - (Siebert and Döll, 2010).However, radiometric measurements provided by remote sensing added to GIS-based climate modelling have proved to be essential for modelling ET because they are the only techniques that allow us to compute it feasibly at both regional (Cristóbal et al., 2005;Kustas and Norman, 1996) and global scales (Mu et al., 2007).Moreover, the use of remote sensing techniques supplements the frequent lack of ground-measured variables and parameters that are required for applying the local models at a regional scale (Sánchez et al., 2008a).
Currently, there are a wide variety of remote sensing models for calculating ET at global or regional scales, such as METRIC (Allen et al., 2007), SEBAL (Bastiaanssen et al., 1998), TSEB (Kustas and Norman, 2000), ALEXI/disALEXI (Anderson et al., 2004), S-SEBI (Roerink et al., 2000), STSEB (Sánchez et al., 2008b) and the B-Method (Jackson et al., 1977;Seguin and Itier, 1983); other methodologies can be found in Schmugge et al. (2002) or Sánchez et al. (2008a).All these theoretical methods used to estimate ET at a regional scale with remote sensing techniques are derived from the energy balance equation based on the principle of energy conservation in a system formed by soil and vegetation.Most of them try to minimize the inputs from ancillary data (often data from ground meteorological stations) in order to make the algorithms more operational at global scales.However, there is currently no agreement on which method is the most appropriate because this often depends on the application purposes.
Most of these methods have been validated in homogeneous covers (crops or natural vegetation) and flat areas, where a single meteorological station record is used to describe the climate conditions of a large area.In these areas, the use of moderate or coarse spatial resolution remote sensing data is enough to obtain accurate daily ET (ET d ) results.However, in more complex and heterogeneous areas, due to the landscape, orographic or climatic variability, a single meteorological record or remote sensing data with coarse spatial resolution may not be accurate enough to calculate the ET d .Operative GIS climate-based techniques can be used at regional scales in both simple and complex areas (Ninyerola et al., 2007;Pons and Ninyerola, 2008) to achieve higher accuracy and provide the input variables for agriculture and natural vegetation evapotranspiration modelling.
The objective of this paper is to evaluate a simple method for computing daily ET using stand-scale sap flow measurements made in Scots Pine (Pinus sylvestris L.) in a heterogeneous Mediterranean mountain area during a three year period.GIS climate-based regional modelling was used instead of a single meteorological measurement to obtain the meteorological inputs (air temperature and solar radiation) in order to evaluate the performance of this technique.Low (TERRA/AQUA MODIS) and moderate (Landsat-TM/ETM+) spatial resolution remote sensing images were used as the remote sensing inputs in order to evaluate their accuracy in a heterogeneous landscape.

Study area
The study plot is located within the Vallcebre research catchments (Gallart et al., 2005;Latron et al., 2010;Llorens et al., 2010), 42 • 12 N, 1 • 49 E, located in the eastern Pyrenees (NE Iberian Peninsula).It has a humid Mediterranean climate, with a marked water deficit in summer.The mean annual temperature at 1260 m is 9.1 • C, and the long term  mean annual precipitation is 862 ± 206 mm, with a mean of 90 rainy days per year.The long term (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) mean annual reference evapotranspiration, calculated using the Hargreaves and Samani (1982) method, was 823 ± 26 mm.Mudstone and limestone substrates are predominant, resulting in clay soils in the first case, and bare rock areas or thin soils in the latter (Gallart et al., 2002).The vegetation in the area is sub-Mediterranean oak forest (Buxo-sempervirentis-Quercetum pubescentis association), but most of the land was terraced and deforested for cultivation in the past, and then progressively abandoned during the second half of the twentieth century.The present landscape is mainly a mosaic of mesophylous grasslands and patches of Scots pine, which colonized old agricultural terraces after they were abandoned (Poyatos et al., 2003).Figure 1 shows the location of the Vallcebre research catchments.
In the case of air temperature and net radiation regionalization, meteorological data from 161 meteorological stations were downloaded from the Catalan Meteorological Service (SMC) web (meteorological data available at http://www.meteocat.com).Figure 1 shows the spatial distribution of these two sources of meteorological data.
A set of 30 TERRA-MODIS images and 27 AQUA-MODIS images and a set of 11 Landsat-7 ETM+ and 10 Landsat-5 TM images from paths 197 and 198, row 31 were selected to perform the ET d modelling of the Scots Pine forest stand from 2003 to 2005.In order to avoid canopy interception and to make sap-flow measurements fully representative of ET canopy we selected clear sky dates where no precipitation was present during at least 15 days before and after the selected day. Figure 2 shows the temporal distribution of the selected dates aggregated by month.
AQUA/TERRA MODIS images were downloaded using the Land Processes Distributed Active Archive Center gateway (https://wist.echo.nasa.gov/∼ wist/api/imswelcome/).We selected three different types of products which contain the remote sensing data to calculate the ET d : MOD11A1 and MYD11A1 (which contain TERRA and AQUA daily land surface temperature, LST, and emissivity respectively), MOD09GHK and MYD09GHK (which contain TERRA and AQUA daily reflectances, respectively), and MOD05 (which contains daily water vapour).Although image time acquisition is different for each satellite, Landsat and TERRA satellites pass over Catalonia at a similar time, between 09:30 and 10:30 LST (local solar time).AQUA passes over the same area, but between 13:00 and 14:00 LST.

The evapotranspiration model
We used the B-method to compute ET d .This methodology is derived from the model proposed by Jackson et al. (1977), which is based on the energy balance equation and has been used or modified for both natural vegetation and crop areas by different authors (Caselles et al., 1992(Caselles et al., , 1998;;García et al., 2007;Sánchez et al., 2007Sánchez et al., , 2008a;;Seguin and Itier, 1983;Vidal and Perrier, 1989).Seguin and Itier (1983) proposed a modified equation that needs net radiation (R n ) and the difference between LST and air temperature at satellite pass (T i ) as input variables: where subindex d is the daily periods, ET and R n are in mm day −1 and both temperatures are in K. Exponent n is a correction for non-neutral static stability that could be assigned to one, as Seguin and Itier (1983) suggested.Since daily integration of soil heat flux is likely to be close to zero, Eq. ( 1) expresses the daily integrated surface sensible heat flux into the atmosphere (Allen et al., 1998;Carlson and Buffum, 1989;Seguin and Itier, 1983).Due to the importance of the B parameter in calculating ET d , we used two approaches to compute B. B can be defined as an exchange coefficient that in Eq. ( 1) represents an average bulk conductance for the dailyintegrated sensible heat flux.This term is related to the sensible heat flux (H ), one of the most difficult variables to determine in the energy balance equation (Bastiaanssen et al., 1998).There are several approaches that use LST directly, such as the parallel resistance model developed by Norman et al. (1995), and the one developed by Caselles et al. (1992), which is adapted for heterogeneous areas and defined by the following equation: where subindex i means instantaneous and (R nd /R ni ) is called the R n ratio.ρC p is the volumetric heat capacity of air (J kg −1 K −1 ) and r * a is the effective aerodynamic resistance.Measurements of effective aerodynamic resistance are not usually easy to obtain; therefore, we considered the aerodynamic resistance of Pinus sylvestris to be equal to 28.1 m s −1 , as determined by Sánchez et al. (2007), because our study area is similar to that of this previous work (J. S. Sánchez, personal communication, 2009).
In addition, B can also be obtained using the simple equation proposed by Carlson et al. (1995), obtained from a soil-vegetation atmosphere transfer model that integrates the main factors on which B depends, such as wind velocity and aerodynamic resistance; therefore, B can also be defined as: where NDVI * is a scaled vegetation index based on the NDVI and is defined as: where subindex p is the image NDVI value for a given pixel, 0 is a bare soil pixel and s is a fully vegetated pixel.

Landsat and TERRA/AQUA image processing
Landsat images were corrected using the methodology proposed by Palà and Pons (1995) that is based on a first-degree polynomial fit that accounts for the relief using a detailed enough Digital Elevation Model (DEM) obtained from the Cartographic Institute of Catalonia (ICC).This correction also requires a set of ground control points (GCP) that were digitized on screen from 2.5 m digital orthophotos (from the ICC).A planimetric accuracy (obtained with an independent set of GCP) of less than 15 m (half pixel) was obtained.Radiometric correction was carried out following the methodology proposed by Pons and Solé-Sugrañes (1994), which allows us to reduce the number of undesired artifacts due to the atmospheric effects or differential illumination that are results of the time of day, the location on the Earth and the relief (zones being more illuminated than others, shadows, etc).The digital numbers were converted to radiances by means of image header parameters, taking into account the considerations presented by Cristóbal et al. (2004) and Chander et al. (2009).
AQUA/TERRA MODIS reflectance and LST and emissivity products were imported, with all the necessary metadata to process them.Before that, images were reprojected to UTM-31 N. The water vapour product was geometrically corrected using HEG-WIN software (http://newsroom.gsfc.nasa.gov/sdptoolkit/HEG/HEGDownload.html).

Air temperature
Different air temperature input variables are needed to compute net radiation LST and ET d : satellite pass air temperature (T i ), daily mean air temperature (T a ) and daily minimum air temperature (T min ).To regionalize air temperature, we applied a multiple regression analysis with spatial interpolation of residual errors of ground meteorological station data using geographical variables as predictors, such as altitude, latitude, or continentality (Cristóbal et al., 2008;Ninyerola et al., 2000Ninyerola et al., , 2007)).Spatial interpolation of the residuals has been computed using the Inverse Distance Weighted interpolation because this interpolator offers better results than other methodologies, at least in the case of air temperature modelling (Ninyerola et al., 2000).Air temperature data were fitted using 60 % of the meteorological ground stations and cross-validated with the remaining 40 %.In these previous works, T i , T a and T min were obtained with an RMSE of 1.8 K, 1.3 K and 2.3 K, respectively.

Land surface temperature (LST) and emissivity (LSE)
In the case of Landsat-5 TM and Landsat-7 ETM+, the LST was calculated with the methodology proposed by Cristóbal et al. (2009), which is based on the radiative transfer equation and needs air temperature and water vapour as input variables, and present a RMSE of about 1 K compared with radiosonde data.The methodology is designed for a wide range of water vapour values (0 to 8 g cm −2 ) to take into account global conditions.The TERRA/AQUA MODIS water vapour product (MOD05) was used as the water vapour source.The air temperature was computed as explained in Sect.4.2.
To compute LSE we used the NDVI Threshold Method proposed by Sobrino and Raissouni (2000) and Sobrino et al. (2008).This methodology uses certain NDVI thresholds to distinguish between bare soil, fully vegetated and mixed pixels.According to the authors it gives an error of 1 % (Sobrino et al., 2008).

Daily net radiation (R nd )
Instantaneous net radiation is usually computed with the energy balance equation as follows: 5) where subindex i means instantaneous, α is the surface albedo, R s↓ is the incoming short wave radiation, σ is the Stephan-Boltzmann constant; ε is the surface emissivity and ε a is the air emissivity.The three terms of Eq. ( 5) regard to incoming net shortwave radiation, incoming longwave radiation and outgoing longwave radiation, respectively.However, B-method needs R nd instead of R ni as input.Therefore, in order to compute R nd we approached the three terms of Eq. ( 5) on a daily basis as follows: where α is the surface albedo, R sd↓ is the daily incoming short wave radiation, L ↓ d is the daily incoming longwave radiation and L ↑ d is the daily outgoing longwave radiation.Albedo (α) was computed using the Liang (2001) methodology in the case of Landsat-5 TM and Landsat-7 ETM+ images, and the method by Liang et al. (1999) in the case of TERRA/AQUA MODIS images.Both methodologies use a weighted sum of visible, near infrared and medium infrared radiation, and according to the authors the error in estimating albedo is less than 2 %.Daily solar radiation (R sd↓ ) was obtained with the methodology proposed by Pons and Ninyerola (2008).Given a digital elevation model, we can calculate the incident solar radiation at each point during a particular day of the year taking into account the position of the Sun, the angles of incidence, the projected shadows, the atmospheric extinction and the distance from the Earth to the Sun at fifteen minute intervals.The diffuse radiation was estimated from the direct radiation and the exoatmospheric direct solar irradiance was estimated with the Page (1986) equation that Baldasano et al. (1994) fitted with information from Catalonia.L ↓ d was computed by means of the methodology proposed by Dilley and O'Brien (1998) that according to the authors shows a RMSE of 5 W m −2 and a R 2 of 0.99 in is computation.
where α, β and γ are 59.38, 113.7 and 96.96, respectively; w is the water vapour, in kg m −2 , T * is 273.16K, w * is 25 kg m −2 and T a is daily mean air temperature.Finally, L ↑ d was modeled by means of the methodology proposed by Lagouarde and Brunet (1993) as follows: where ε is the land surface emissivity and R is defined as: where σ is Stefan-Boltzmann constant (5.67 × 10 −8 W K −4 m −2 d −1 ), T is the difference between LST and T a at satellite pass, t, (both in K), T min (K), α = 1.13,D is the time difference between sunset and sunrise; and τ = 24 (for a 24 h period).

B parameter
As we explained in Sect.3, B parameter was calculated with two approaches: the R n ratio and NDVI.In the R n ratio approach, we used two ways to compute the parameter: (1) a regional R n ratio (hereafter referred to as the B-R n ratio regional) with data from 13 meteorological stations of the SMC meteorological network, and (2) a local R n ratio (hereafter referred to as the B-R n ratio local) with data from the meteorological station above the Scots pine stand in the Vallcebre research catchments.We used these two data sources to evaluate whether a regional measurement of the R n ratio provides similar results as a local measurement.
In the NDVI approach (hereafter referred to as the B-NDVI), Carlson et al. (1995) suggested selecting NDVI values depending on the study area.In our case, bare soil and fully vegetated NDVI values were set to 0.1 and 0.7 for the entire dataset, as these values were realistic enough to simulate bare soil and full vegetation conditions over the study area.

Sap flow measurements and upscaling to stand transpiration
We compared remote sensing daily evapotranspiration estimates with sap flow measurements upscaled to stand transpiration.Sap flow density in the outer xylem was measured with 20 mm long heat dissipation probes constructed according to Granier (1985); 15-min averages of data collected every 10 seconds were stored in a datalogger (DT 500, DataTaker, Australia).Heat dissipation gauges were installed at breast height on the north-facing side of 12 Scots pine trees and were covered with reflective insulation to avoid the influence of natural temperature gradients in the trunk.Sap flow density measured in the outer xylem was corrected for radial variability in sap flow using correction coefficients derived from radial patterns of sap flow within the xylem measured with a multi-point heat field deformation sensor (Nadezhdina et al., 2002).A gravimetric analysis of wood cores was carried out to estimate sapwood depths in a sample of Scots pine trees, and a linear regression was obtained between the basal area and sapwood area of individual trees.Stand transpiration was then calculated by multiplying the average sap flow density within a diametric class by the total sapwood area of trees in that class.Instantaneous values (15-min averages) were then summed to compute daily stand transpiration.Further details on the methodology and results of sap flow measurements used in this study can be found in Poyatos et al. (2005Poyatos et al. ( , 2008)).

B parameter
The B parameter showed different behaviour depending on the approach used.The B-R n ratio local had a mean and standard deviation (s.d.) of 6.9 and 3.2 W m −2 respectively, in the case of Landsat dates (see Fig. 2), and a mean and s.d. of 10.8 and 2.2 W m −2 respectively, in the case of TERRA/AQUA dates.The B-R n ratio regional displayed a mean and s.d. of 9.5 and 2.3 W m −2 respectively, in the case of Landsat dates, and a mean and s.d. of 12.8 and 2.4 W m −2 respectively, in the case of TERRA/AQUA dates.Finally, B-NDVI showed a mean and s.d. of 11.9 and 2.6 W m −2 respectively, in the case of Landsat dates, and a mean and s.d. of 12.6 and 2.6 W m −2 respectively, in the case of TERRA/AQUA dates.B-NDVI was similar in Landsat and TERRA/AQUA dates, but not in the B approach using the R n ratio, especially in the case of the B-R n ratio local.While on winter and autumn dates the B-R n ratio local had small values (positive or negative) close to 0, B-NDVI tended to show higher positive values.For example, B computed on 11 January 2005, using the local R n ratio gave a negative value close to 0 W m −2 , whereas in the case of B-NDVI it was 11.8 W m −2 .During these seasons we would expect low B values due to the energy budget; therefore, this suggests that B-NDVI could be less sensitive in winter and autumn situations than the B-R n ratio.
In the case of the B-R n ratio, the R n ratio is usually obtained from a net radiation sensor over the study area.Some authors have used a constant value of 0.3 ± 0.02 (Seguin and Itier, 1983;Kustas et al., 1990;García et al., 2007) because most of the dates used in these studies were in spring or summer and over crop areas.However, we found that our local (Vallcebre research catchments) and regional (SMC meteorological stations) R n ratios varied over the year (see Fig. 3).The local and regional R n ratios for the Landsat/TERRA satellite pass had an annual mean (from 2003 to 2005 period) of 0.16 ± 0.05 (mean and s.d.) and 0.22 ± 0.03 respectively, and in the case of the AQUA satellite pass, an annual mean of 0.17 ± 0.05 and 0.21 ± 0.02 respectively.In addition, the R n ratio varied little from 09:00 to 14:00 in our study area, and thus was useful in Landsat and TERRA/AQUA ET d modelling.Therefore, we used a daily R n ratio instead of a constant R n ratio.This is in agreement with other authors who also reported a similar R n ratio behaviour (Sánchez et al., 2007;Sobrino et al., 2005;Wassenaar et al., 2002).R n ratio values reported in these studies are similar to the regional R n ratio computed in our study area because our value was obtained at meteorological stations at a similar altitude as those in the literature.However, the local R n ratio values are lower, which shows that this ratio does not only vary with latitude as Sánchez et al. (2007) suggests, but also with altitude.
Hydrol.Earth Syst.Sci., 15, 1563Sci., 15, -1575Sci., 15, , 2011 www.hydrol-earth-syst-sci.net/15/1563/2011/Further research into R n ratio modelling in mountainous areas is therefore needed.Moreover, it is worth remarking that the local R n ratio obtained in the Scots pine stand displayed a different pattern to the regional R n ratio because it had negative values on winter days when the R n budget is negative, which often occurs in mountainous areas (Barry, 2001).ET d models do not usually predict this situation because they are generally applied in spring and summer and in relatively flat and low altitude areas.

Daily net radiation validation
The results show a good agreement between the R nd measured in the Scots pine stand and the R nd obtained using the Landsat regional model with an R 2 of 0.89 and an RMSE of 22 W m −2 (see Table 1).The R nd derived from TERRA/AQUA MODIS showed a similar RMSE but lower R 2 (0.77 and 0.73 respectively; Table 2).It is worth noting that the proposed R nd model developed with regional variables, such as R sd↓ , LST and T a , makes it possible to approximate this variable over large areas with a high level of accuracy.

Landsat
In the ET d validation, we obtained a test R 2 of 0.84 when the B parameter was estimated using a regional R n ratio computed with the SMC meteorological stations (B-R n ratio regional), 0.84 when the B parameter was estimated using a local R n ratio computed with the Scots pine stand meteorological station (B-R n ratio local), and 0.82 when the B parameter was estimated using the NDVI approach, B-NDVI.It is interesting to note that for the ET d models used, the minimum values were always negative.This mainly happens on winter dates when the R n ratio is also negative.Therefore, on winter dates this methodology should only be used on days when the R n budget is positive.Errors close to 1 mm day −1 were obtained for the RMSE.Taking into account the range of ET d values observed in the studied Scots pine stand (from 0.5 to 2.7 mm day −1 ), we cannot conclude that the model provides optimal results.When the R nd ratio is negative during winter, the ET d yields negative values and the model does not perform well.Again, it is worth noting that ET d models are usually validated on spring or summer dates (Chiesi et al., 2002;Nagler et al., 2005Nagler et al., , 2007;;Sánchez et al., 2007Sánchez et al., , 2008a;;Verstraeten et al., 2005;Wu et al., 2006) when the daily R n budget is positive.Our attempt to also estimate ET d during autumn and winter has shown the limitations of the method and how ET d modelling needs to be further improved, especially in forest areas.
We obtained a better mean RMSE for the different models when only those dates with a positive R nd ratio were selected, which ranged from 0.5 to 0.7 mm day −1 with a similar R 2 (see Tables 2-3 and Fig. 4).Of the different approaches used to compute the B parameter, the best results were obtained using the local R n ratio and the NDVI approaches, with a RMSE of 0.5 and 0.6 mm day −1 , respectively, and an estimation error of about ±30 %.Indeed, the regional R n ratio yielded a higher RMSE and estimation error  3 B-Rn ratio regional is the B approach using the regional Rn ratio, B-Rn ratio local is the 4 B approach using the local Rn ratio and B-NDVI is the B approach using the NDVI. 5 Fig. 4. Relationship between daily evapotranspiration (ET d ) calculated from sap flow measurements modelled using Landsat data and the B parameter approach in mm day −1 .B-R n ratio regional is the B approach using the regional R n ratio, B-R n ratio local is the B approach using the local R n ratio and B-NDVI is the B approach using the NDVI. of about 0.7 mm day −1 and ±38 % respectively; this could be explained by the differences in the R n ratio estimation.The regional R n ratio was computed from the data from the SMC meteorological stations, which are designed for crop assessment and are located in areas at low to medium heights (from 0 to 500 m).Our study area is located at 1250 m; therefore, the regional R n ratio conditions are not representative of our study area.However, it is important to note that optimal R n ratio values are difficult to obtain because it would be necessary to have a meteorological network with R n instruments distributed at different altitudes in diverse landscapes.Moreover, R n instruments are usually found in agrometeorological networks but not very often over forest areas.Although the two B parameter approaches (B-R n ratio local and B-NDVI) obtained similar results, the main advantage of the NDVI approach is easily implemented, when realistic values of NDVI thresholds to compute NDVI * are selected, to compute ET than B-R n ratio local or regional because that require intensive local calibration.So, if a well-balanced regional R n ratio is not available due to limitations in the meteorological networks, the NDVI approach is preferable for computing the B parameter at regional scales in an operative way.In all cases, the models tended to overestimate ET d , showing higher values in the case of the regional R n ratio and lower values in the case of the NDVI approach.
In this study, we are strictly comparing evapotranspiration with stand transpiration of the dominant tree species.As the understory in the studied stand is very poor, the only other contribution comes from soil evaporation, with typical rates of 0.1 to 0.5 mm day −1 (Poyatos et al., 2007).These values are consistent with the systematic bias between sap flowderived transpiration and the ET d models (see Fig. 4).
In addition, it should be stressed that the difficulty of obtaining the effective aerodynamic resistance and the use of a constant value for the analyzed period may have introduced more variability into our analysis, thus increased the error in the ET d models.

TERRA/AQUA MODIS
TERRA and AQUA ET d validation did not obtain the same results as Landsat (see Tables 2-3 and Figs.5-6).In both cases, ET d validation showed a higher RMSE, between 1.8 and 2.4 mm day −1 , a low R 2 , between 0.03 and 0.07 and estimation error of about ±57 and 50 %, respectively.Despite air temperature models also showing good validation results, one possible source of error is the LST.Therefore, it seems that although TERRA/AQUA MODIS LST provides good results in R nd modelling, this is not the case in ET d modelling.
The study area does not cover a 1000 m × 1000 m pixel, and therefore the remote sensing data, especially the LST, is less representative than a Landsat pixel of 60 (ETM+) or 120 m (TM).As we mentioned in Sect.2.1, the study area has high spatial heterogeneity, mosaic of afforestation patches overgrowing old agricultural terraces (Poyatos et al., 2003), at smaller scales than the coarse TERRA/AQUA MODIS pixel, which makes it difficult to validate the ET d model results.However, it has to be noted that validating a TERRA/AQUA MODIS in a heterogeneous mountain landscape is not easy due to the extensive instrumentation needed to measure the energy flux in each of the landscape covers.Therefore, it seems that the use of TERRA/AQUA MODIS images themselves in this type of landscape is not enough in mm•day -1 .B-R n ratio regional is the B approach using the regional R n ratio, B-R n ratio 5 local is the B approach using the local R n ratio and B-NDVI is the B approach using the 6 NDVI.
7 Fig. 5. Relationship between daily evapotransporation (ET d ) calculated from sap flow measurements and modelled using TERRA MODIS data and the B parameter estimation in mm•day −1 .B-R n ratio regional is the B approach using the regional R n ratio, B-R n ratio local is the B approach using the local R n ratio and B-NDVI is the B approach using the NDVI.
to accurately map ET on a daily basis.In order to improve the ET d results using coarse resolution images, downscaling techniques such as those in Anderson et al. (2004) are required.
There are very few studies in the literature that monitor ET d at both high spatial and temporal resolutions during an annual period in a forest area, especially using a large number of Landsat images.In addition, most of the studies to date have dealt with environments subjected to in mm•day -1 .B-R n ratio regional is the B approach using the regional R n ratio, B-R n ratio local is the B approach using the local R n ratio and B-NDVI is the B approach using the NDVI.
Fig. 6.Relationship between daily evapotransporation (ET d ) calculated from sap flow measurements and modelled using AQUA MODIS data and the B parameter estimation in mm day −1 .B-R n ratio regional is the B approach using the regional R ratio, B-R n ratio local is the B approach using the local R n ratio and B-NDVI is the B approach using the NDVI.
only mild water stress, such as riparian forests, crops or boreal stands.For example, Wu et al. (2006) reported an RMSE of 0.6 mm day −1 in a tropical forest using only one Landsat image.They compared their results with estimates from the literature due to the difficulty in validating these kinds of regions using sap flow measurements.Nagler et al. (2005Nagler et al. ( , 2007) )  Overall, for Landsat ET d modelling, our results are in agreement with the studies in the literature, as we obtained an uncertainty of about 30 %.One positive point of the results of the Landsat ET d models across different seasons is the robust ET d estimation under varied conditions of water availability, as the studied stand undergoes different degrees of water stress in spring, summer and autumn (Poyatos et al., 2008).However, in the case of MODIS ET d modelling, the validation shows a higher RMSE, which suggests that higher spatial resolution is needed for heterogeneous areas.
In addition, it is worth noting that implementing regional models for calculating T a , LST and R sd↓ , as inputs in both R n and ET d modelling has provided good results and made it possible to compute these variables at regional scales with similar accuracy to that in the literature.
It is important to note, however, that we have found some limitations in ET d modelling in a mountainous forested area that should be addressed in the future in order to monitor this variable in an operational way.Further work to improve the described methodology should include: (i) the validation of a multi-scale remote sensing model (Anderson et al., 2004) for disaggregating regional fluxes to micrometeorological scales.This would allow ET d to be monitored on a daily basis instead of on a 16-day basis thanks to the TERRA/AQUA temporal resolution; (ii) the implementation of methodologies for calculating aerodynamic resistance, such as those described by Norman et al. (1995) and Sánchez et al. (2008a,b).

Conclusions
The B-method has been used to estimate daily evapotranspiration (ET d ) in a Scots pine stand in a mountainous Mediterranean area, obtaining an estimation error of ±30 % (corresponding to 0.5-0.7 mm day −1 ) using medium spatial resolution imagery, Landsat-5 TM and Landsat-7 ETM+, and the different approaches presented.These results are in agreement with recent studies that used a similar spatial resolution.However, when lower spatial resolution was used (TERRA/AQUA MODIS) the results showed larger errors, 1.9 and 1.7 mm day −1 respectively.
The R n ratio emerged as an important parameter to be considered when the B-method is used.Although this ratio is close to 0.3 in spring and summer months, this value is not appropriate for winter and autumn because when the R n ratio is negative (negative R n budget) the B-method does not provide a realistic ET d .Further research is therefore needed to estimate this parameter in these conditions.
The best ET d results were obtained using a local R n ratio approach to calculate the B parameter, followed by the method using NDVI.The regional R n ratio resulted in larger errors, which means that if a well balanced meteorological network (with R n sensors) is not available, the NDVI approach is preferable for calculating the B parameter at a regional scale in an operative way.
Regional input variables for calculating ET d , such as R sd↓ , LST and T a , performed well, making possible to compute it at a regional scale with a good level of accuracy.
Hydrol.Earth Syst.Sci., 15, 1563Sci., 15, -1575Sci., 15, , 2011 www.hydrol-earth-syst-sci.net/15/1563/2011/ Finally, using a large number of remote sensing images that are well distributed over the analyzed period, especially in the case of Landsat, allowed us to better understand the limitations of the methodologies and how to address the further improvement of ET d modelling, especially in forest areas.

Fig. 1 .
Fig. 1.Location of SMC meteorological stations and Vallcebre research catchments in Universal Transversal Mercator (UTM) projection (UTM coordinates are expressed in km).The white dots are meteorological stations from the SMC that include air temperature sensors, the black dots are meteorological stations from the SMC that include net radiation sensors, and the black triangle indicates the Vallcebre research catchments.Panel (A) is the Landsat-TM LST of 1 July 2003 and panel (B) is the TERRA MODIS LST of 10 July 2003 of the Vallcebre research catchments (black triangle).The red square represents a Landsat-TM thermal band pixel (120 m) and the yellow square represents a TERRA MODIS thermal band pixel (1000 m).In panels (A) and (B) the white dot is the Scots pine stand.

Figure 2 .
Figure 2. Monthly temporal distribution of satellite data (clear sky and without bow tie effect: an artefact of the arrangement of sensors on the MODIS instrument, in which the scans are partially overlapping at off nadir angles) used during the period 2003-2005.

Fig. 2 .
Fig. 2. Monthly temporal distribution of satellite data (clear sky and without bow tie effect: an artefact of the arrangement of sensors on the MODIS instrument, in which the scans are partially overlapping at off nadir angles) used during the period 2003-2005.

Figure 3 . 4 Fig. 3 .
Figure 3. Mean monthly local and regional R n at the Landsat/TERRA and AQUA 2 1

Figure 5 .
Figure 5. Relationship between daily evapotransporation (ET d ) calculated from sap flow 3

Figure 6 .
Figure 6.Relationship between daily evapotransporation (ET d ) calculated from sap flow measurements and modelled using AQUA MODIS data and the B parameter estimation modelled ET d using 8 Landsat-5 TM and about 90 MODIS dates in a cottonwood plantation in riparian corridors of the western US during the July-August period in 2005, and obtained an uncertainty in modelled ET d of 20-30 %.Verstraeten et al. (2005) also reported an uncertainty of about 27 % in instantaneous ET modelled with NOAAimagery and validated using EUROFLUX data during the growing seasons of European forests, from March to October.Sánchez et al. (2007) modelled ET d using MODIS in a homogeneous Pinus sylvestris stand in the boreal region, and reported an RMSE of 0.81 mm day −1 and an uncertainty of about 30 % in ET d compared to eddy-covariance data.With a more demanding method in terms of ancillary data needs, FOREST-BGC, Chiesi et al. (2002) reported a mean RMSE of 0.4 mm day −1 introducing LAI derived from 10-day composites of NOAA images in two oak stands.More recently, Sánchez et al. (2008a) obtained a value of 0.7 mm day −1 in different coniferous, broad-leaf and mixed forests in the Basilicata region with three Landsat-5 TM and ETM+ images from spring and summer.

Table 1 .
(a): Descriptive statistics of daily net radiation (Rn d ) and daily evapotranspiration (ET d ) measured over the Scots pine stand, and Landsat Rn and ET d modelled using 3 methods: B-R n ratio regional, B-R n ratio local and B-NDVI (see text).(b) Model validation results.s.d. is standard deviation, RMSE is root mean square error and MBE is mean bias error.

Table 2 .
(a): Descriptive statistics of daily net radiation (Rn d ) and daily evapotranspiration (ET d ) measured over the Scots pine stand, and TERRA MODIS Rn and ET d modelled using 3 methods B-R n ratio regional, B-R n ratio local and B-NDVI (see text).(b) Model validation results.s.d. is standard deviation, RMSE is root mean square error and MBE is mean bias error.

Table 3 .
(a): Descriptive statistics of daily net radiation (Rn d ) and daily evapotranspiration (ET d ) measured over the Scots pine stand, and AQUA MODIS Rn and ET d modelled using 3 methods: B-R n ratio regional, B-R n ratio local and B-NDVI (see text).(b) Model validation results.s.d. is standard deviation, RMSE is root mean square error and MBE is mean bias error.