Two-source energy balance modeling of evapotranspiration in Alpine grasslands

Remote Sensing of Environment 209 (2018) 327–342

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Review

Two-source energy balance modeling of evapotranspiration in Alpine grasslands

T

M. Castellia,*, M.C. Andersond, Y. Yangd, G. Wohlfahrtc, G. Bertoldib, G. Niedristb, A. Hammerlec, P. Zhaoc, M. Zebischa, C. Notarnicolaa a b c d

Institute for Earth Observation, Eurac Research, via Druso 1, Bolzano 39100, Italy Institute for Alpine Environment, Eurac Research, via Druso 1, Bolzano 39100, Italy Institute for Ecology, University of Innsbruck, Innsbruck, Austria U.S. Department of Agriculture, Agricultural Research Service, Hydrology and Remote Sensing Laboratory, Beltsville, MD, USA

A B S T R A C T This work aims to assess a diagnostic approach which links evapotranspiration (ET) to land surface temperature (LST) measured by thermal remote sensing in the Alps. We estimated gridded ET, from field (30 m) to regional (1 km) scales, and we performed a specific study on grassland ecosystems in the Alps in South Tyrol (Italy), to evaluate the model sensitivity to different types of land management. The energy balance model TSEB ALEXI (Two Source Energy Balance Atmosphere Land EXchange Inverse) was first applied to Meteosat satellite data. Then ET was estimated by the flux disaggregation procedure DisALEXI, driven by MODIS and Landsat LST retrievals, which has never been applied before in a mountain region. We evaluated the model against eddycovariance (EC) measurements from established stations in the Alps, and analyzed the main limitations which affect the model performance in mountainous regions. The TSEB model, applied in plot-scale mode using towerbased meteorological and LST input data, performed well with errors in daytime (6–18 UTC+1) latent heat flux around 30–60 W m−2 in comparison with flux measurements corrected for the lack of closure in the energy balance. For landscape ET retrievals, while Landsat resolution (30 m) is preferable for capturing small-scale heterogeneity in landscape moisture conditions, and for direct comparison with tower fluxes, persistent cloud cover resulted in no clear Landsat scenes during the study period. MODIS-based retrievals at 1 km resolution are too coarse to resolve the flux tower footprint in this complex landscape, yielding discrepancies of 100 W m−2 in model-measurement comparisons. Still, MODIS DisALEXI partitioning of the energy budget was reasonable and enabled to detect evaporative stress at regional scale expressed as the ratio between actual and potential ET, fPET. We evaluated fPET in comparison with a crop stress index based on cumulative air temperature and precipitation at different stations in the study area, and investigated ability to capture differential responses between managed and unmanaged grasslands. Results show that in the Alps i) moderate resolution thermal data can be used to monitor evaporative stress at the regional scale; ii) the spatial-temporal evolution of ET can be characterized from MODIS and Landsat thermal data with limitations which are due to the low availability of clear-sky scenes and to the small-scale (∼10 m) changes in soil moisture, topography and canopy density, which control ET patterns in mountainous regions; iii) solar radiation and leaf area index are critical variables which strongly affect the accuracy of the modeled energy fluxes.

1. Introduction Among land–atmosphere fluxes, evapotranspiration (ET) plays a key role in the global biogeochemical cycling since it links the water, energy and carbon cycles. This makes ET an important indicator for monitoring drought conditions. In the European Alps, summer droughts are getting increasingly frequent with climate change mainly due to a

*

shift of precipitation from summer to winter and higher ET in summer due to higher temperatures (Haslinger et al., 2015). Dry summer conditions are already a frequent phenomenon in inner-alpine dry valleys such as our case studies Val Venosta/Vinschgau and Val Mazia/Matschertal. Furthermore, mountainous regions are especially vulnerable because temperature increases are expected to be faster than the global

Corresponding author at: via Druso 1, Bolzano 39100, Italy. E-mail address: [email protected] (M. Castelli).

https://doi.org/10.1016/j.rse.2018.02.062 Received 29 July 2016; Received in revised form 30 November 2017; Accepted 22 February 2018 0034-4257/ © 2018 Elsevier Inc. All rights reserved.

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

The analysis presented in this paper is mainly concentrated on grassland, because it is part of a project which encompasses both microwave and thermal infrared remote sensing. In parallel with ET estimation, soil moisture is modeled from radar data (Pasolli et al., 2015b) by machine learning algorithms. Due to the reduced penetration capacity of microwaves, soil moisture is estimated only on bare soil and areas with limited vegetation cover. Ground measurement stations, including eddy covariance systems and sensors measuring soil moisture, have been established in grassland areas. The structure of this paper is as follows: Section 2 describes the study area in which validation instrumentation and ground-based measurements are located, Section 3 summarizes the characteristics of the TSEB ALEXI/DisALEXI model, the model input used for the simulations at the plot and at the regional scale, and the processing of EC data, Section 4 illustrates the main results of this study, including the evaluation of the plot scale and of the remote sensing-based simulations, a sensitivity analysis on TSEB and the derivation and assessment of an evaporative stress index, and in last section the conclusions derived from the outcome of the present analysis are summarized and discussed, and possible continuations of the present study are proposed.

land average, and decrease in snow fall may have additional impacts on the hydrological cycle (Gobiet et al., 2014). The uncertainty on future climate requires new ways of managing agriculture that are adaptable to the occurrence of extreme events. In this context of climate transitions, monitoring ET is important to assist policy makers as they develop strategies to reduce the vulnerability of agriculture and forestry to climate change. At the watershed scale, ET represents the most important flux of the water balance next to precipitation, but it is the most challenging variable to model due to the heterogeneity of soil and vegetation characteristics. In addition, it is very difficult to accurately simulate the stomatal conductance, which is a key parameter for modeling ET. These problems are exacerbated in mountainous regions, due to the small scale variability of land-cover, soil type, and incoming solar radiation. Various approaches have been developed for modeling ET, which can be classified as diagnostic or prognostic models. Prognostic models are based on the water budget and need spatially distributed precipitation in input. In such models, ET is regulated by the available water in the soil, the radiation input, climatic variables, soil and vegetation properties. In contrast, diagnostic models link ET to Land Surface Temperature (LST) derived from thermal infrared (TIR) remote sensing. The main advantage of the latter approach is that it does not require spatially distributed measurements of precipitation and of the soil hydraulic properties. This is a huge benefit in regions characterized by data scarcity and by strong heterogeneity, such as mountainous areas in the Alps. Thermal remote sensing data available in the Alps include Meteosat, MODIS and Landsat. Meteosat has a nominal resolution of 3 km and transmits images every 15 min, MODIS LST has a resolution of 1 km and is available daily, and Landsat has a resolution of 100–120 m and has a revisit time of 16 days. Clearly, each dataset has advantages – Meteosat and MODIS offering a high temporal resolution and Landsat a high spatial resolution – and an integrated approach has been developed for an optimized exploitation of these data for mapping daily ET at the resolution of Landsat (Cammalleri et al., 2014a). Different methods have been proposed for estimating ET from LST derived from TIR imagery. Due to the uncertainties in the LST retrieval process, such methods adopt strategies for limiting the resulting error. Some of them exploit the spatial variability of LST. It is the case of SEBAL (Surface Energy Balance Algorithm for Land, Bastiaanssen et al., 1998) and METRIC (Mapping Evapotranspiration with Internalized Calibration, Allen et al., 2007), which scale the thermal image after having identified the extreme pixels, corresponding to limiting and nonlimiting water availability in the soil. Some other models, in contrast, like ALEXI (Atmosphere–Land Exchange Inverse) (Anderson et al., 1997, 2007), adopt a time-differencing scheme (Kustas et al., 2001) which is based on the temporal variability of LST, and not on its absolute value. Specifically, the morning LST rise is measured from two images collected by a geostationary satellite during the morning. The advantage with respect to models like SEBAL and METRIC is that the flux estimation over large regions can be automated, and the spatial variability of the meteorological conditions over the modeling domain is considered. SEBAL and METRIC, in fact, calibrate the energy balance at two extreme conditions (dry and wet) relying on the availability of high quality local weather data and on the operator's ability to select appropriate cold and hot pixels. The model TSEB DisALEXI has been already exploited in many different regions all over the world. However, it has been not optimized yet for mountainous areas, where vegetation characteristics (LAI, fractional cover, height of the canopy) change rapidly, agricultural areas are irrigated in a fragmented way, and topography is complex. The aim of this work is to assess the model performances and identify the most critical points for understanding to which extent can it be applied in its status to mountainous regions, and where a major effort is required for adapting the model to areas characterized by complexity in topography and heterogeneity in land-cover.

2. The study area The area of interest is located in the central Italian Alps (Fig. 1) and includes the Provinces of Trento and Bolzano/Bozen. Precipitation data collected at different meteorological stations in both the Provinces (Fig. 1) were used for evaluating the thermal stress index derived from MODIS-based ET estimates. Model validation results are based on the data collected in the emerging Long-Term Ecological Research (LTER) site Val Mazia/ Matschertal (Bertoldi et al., 2014), in the North West of the Province of Bolzano/Bozen, with an area of approximately 100 km2. The predominant land use types are meadow, pasture, and forest. Due to its inner-alpine location, the area is relatively dry (525 mm at 1600 m a.s.l. mean of the years 1920–2010), hence meadows are additionally irrigated up to 6 times per year. Mean annual temperature (on the meadow) is 6.7 °C, and the mean snow cover duration is about 82 days per year. Two micrometeorological measurement stations have been set up in a distance of 500 m: one in an irrigated meadow and one in a rainfed pasture. Both areas are south-west exposed and slightly inclined. Vegetation in the meadow has been characterized as Trisetum flavescens and is mown twice a year (end of June, begin of September). Pasture vegetation belongs to the Festuca valesiaca and is lightly grazed by cattle from mid-June to mid-October. 3. Data and method 3.1. Micrometeorological measurements The two grassland measurement stations in Val Mazia/Matschertal are equipped with a 4 component net radiometer (Hukseflux and Kipp & Zonen CNR-1), soil heat flux plates at 7 cm depth (Hukseflux), soil temperature at 5 and 10 cm depth (Campbell Scientific TCAV) and further standard-meteorological sensors for measuring air temperature and humidity (HMP45C, Campbell Scientific, USA), wind (Windsonic, Gill, UK) and precipitation (Young). Horizontally measured net radiation was corrected for slope and aspect of the underlying surface as detailed in Wohlfahrt et al. (2016). The net ecosystem exchange of sensible and latent heat was quantified using the eddy covariance (EC) method (Aubinet et al., 2000; Baldocchi et al., 1988) at the meadow and pasture using identical instrumentation and post-processing methods. Measurements at the meadow included the period 2012–2014, at the pasture 2014–2015. The three wind components and the speed of sound were measured with a three-dimensional sonic anemometer (CSAT3, Campbell 328

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

Fig. 1. a) Provinces of Bolzano/Bozen and Trento, Italy (Wikimedia Commons). b) Position of measurement stations in the Provinces of Bolzano/Bozen (North) and Trento (South). c) MODIS grid overlapped to the Google maps image of the study area in Val Mazia/ Matschertal. The triangle represents the position of the measurement station in the meadow, and the circle represents the one of the station in the pasture.

is not generally worse for mountain grasslands in complex topography (e.g. Hammerle et al., 2007; Hiller et al., 2008). Considering that when the energy imbalance of observed fluxes is too high, the utility of LE and H for model validation or calibration is greatly reduced, we excluded measurements from our analyses when absolute closure (RN-G-LE-H) exceeded 300 W m−2, and when relative closure ((RN-G-LE-H)/RN) exceeded 0.4. We also tested a stricter criterion, i.e. we excluded measurements when the discrepancy in the energy balance closure, D = (H + LE) / (RN-G), was below 0.7, but this did not significantly affect the results in terms of model performance, and thus we stayed with the above criteria. In addition, for confronting, for calibration or validation, measured components of the surface energy balance with simulations based on a model (TSEB in this case) which assumes the energy balance to be closed, closure needs to be forced with the measured data (Wohlfahrt et al., 2009). To this end, we explored two commonly used options (Twine et al., 2000; Wohlfahrt et al., 2009; Wohlfahrt and Widmoser, 2013): first, we assigned the residual energy to the measured latent heat flux only, referred to as Res in the following; second, we assumed that the Bowen ratio, the ratio between the sensible to latent heat flux, was measured correctly by the eddy covariance method, and distributed the residual energy in proportion to the Bowen ratio to the sensible and latent heat flux (referred to as BR). In case of attributing all the residual of the energy budget to LE, which represents a sort of “worst-case scenario”, the corrected LE can be expressed as follows:

Scientific, USA) mounted 1.7 (meadow)/2 (pasture) m vertically above ground. Water vapor mole densities were measured with an open-path infrared gas analyzer (Li-7500, LiCor, USA), which was displaced by 0.1 m laterally and 0.2 m longitudinally compared to the sonic anemometer in order to minimize flux loss due to sensor separation and to avoid disturbance of the wind field. A data logger (CR3000, Campbell Scientific, USA) acquired the serial data from both instruments at 20 Hz. Post-processing of the raw data was conducted using the free software EdiRe (University of Edinburgh, UK). Half-hourly average fluxes of latent and sensible heat were calculated as the covariance between the turbulent fluctuations of the vertical wind speed and the water vapor density and sonic temperature, respectively, after subtracting the time series arithmetic average and applying a planar fit rotation (Wilczak et al., 2001) to the wind data. Frequency response corrections accounting for both high- and low-pass filtering were applied following Moore (1986) and Aubinet et al. (2000) using sitespecific model cospectra as in Wohlfahrt et al. (2005). The buoyancy flux was converted to the sensible heat flux according to Schotanus et al. (1983) and latent heat fluxes were corrected for the effects of density fluctuations based on Webb et al. (1980). No corrections were applied for the self-heating of the open-path infrared gas analyzer, as these are known to be negligible during the warm season (Burba et al., 2008; Haslwanter et al., 2009). The net exchange of sensible and latent heat was then calculated as the sum of the respective corrected eddy covariance and the storage flux, which was calculated as the vertically integrated rate of change of latent/sensible heat at the measurement height. Generally the combination of EC measurements of latent heat (LE) and sensible heat (H) with independent measurements of net radiation (RN) and soil heat flux (G) generates a lack of closure of the energy balance for several reasons discussed e.g. in Twine et al. (2000) and in Wohlfahrt and Widmoser (2013). A detailed analysis of the causes of the energy imbalance is beyond the scope of the present study and in fact is an outstanding micrometeorological problem, which is discussed extensively in the community (see review by Foken, 2008). Despite it is clear that the assumptions underlying the eddy covariance method (Aubinet et al., 2000) are more likely to be violated in complex mountainous terrain as compared to horizontally homogeneous flat terrain, here we note that the lack of energy balance closure in this study falls well into the range observed across a wide range of FLUXNET sites (Stoy et al., 2013), and that the energy balance closure

Res LEcorr = RNobs − Gobs − Hobs .

(1)

On the other hand, if we assume that the EC correctly measures the BR BR / LEcorr Bowen ratio (Hobs/LEobs), i.e. Hobs / LEobs = Hcorr , and if A is the available energy, defined as A = RNobs − Gobs, and ϵ is the unclosure of the energy budget, defined as ϵ = A − (LEobs + Hobs), flux corrections and corrected fluxes can be expressed as follows:

dH = [BR (LEobs + Hobs ) − Hobs]/(1 + BR)

(2)

dLE = ϵ − dH

(3)

BR LEcorr = LEobs + dLE

(4)

BR Hcorr = Hobs + dH

(5)

where the subscripts obs and corr refer to the uncorrected and corrected 329

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

fluxes, and dLE and dH are the corrections applied to H and LE. In the following, we refer to the uncorrected, i.e. as measured, fluxes as NC (no closure).

where cG varies during the day and according to the soil type and soil moisture conditions. However, cG can be assumed constant during the midmorning-midday hours (Santanello and Friedl, 2003). The TSEB model described so far has been implemented in the algorithm ALEXI (Anderson et al., 1997, 2007), developed for exploiting operationally thermal remote sensing data. To compensate for errors in air temperature and in LST retrieval, ALEXI runs the TSEB model two times, soon after sunrise and before local noon, during the development of the ABL. Land–atmosphere interactions (effect of TA rise on H) are simulated by an ABL model, combined with the TSEB model, which couples the rise in air temperature above the canopy and the growth of the ABL with the sensible heat which is transferred from the surface to the atmosphere during the time interval between the two thermal measurements. Thus, the ABL model internally calculates air temperature at the mixing height, which is not required as input. Thanks to the coupling of the TSEB with the ABL model, ALEXI elaborates only time-differential temperature signals, thus reducing the model sensitivity to errors in the estimation of LST (Kustas et al., 2001). The flux disaggregation scheme DisALEXI was applied to generate consistent flux maps from the MSG (3 km) scale to the MODIS (1 km) and Landsat (30 m) scale. In DisALEXI, air temperature diagnosed by ALEXI is used as upper boundary condition for running the TSEB on single maps of LST and LAI derived from MODIS and Landsat. Daytime total H of ALEXI is used to guarantee the consistency between ALEXI and DisALEXI fluxes. Specifically TA for DisALEXI pixels is iteratively modified until the average daytime H of the DisALEXI pixels which cover the ALEXI pixel matches the one of ALEXI (Anderson et al., 2012).

3.2. Two-source energy balance model description In this work, we exploit the modeling framework ALEXI, which couples a two-source land-surface energy balance model (TSEB) with an atmospheric boundary layer (ABL) model. The TSEB partitions the observed directional radiometric temperature, TRAD(θ), into soil and canopy temperature (TS and TC) based on the vegetation cover fraction at the thermal sensor view angle, f(θ): 1

TRAD (θ) = {f (θ) TC4 + [1 − f (θ)] TS4 } 4 .

(6)

For a homogeneous canopy with a spherical LAI distribution, f(θ) can be approximated as follows:

f (θ) = 1 −

exp ⎛ ⎝ ⎜

−0.5Ω(θ) LAI ⎞ cos(θ) ⎠ ⎟

(7)

where Ω(θ) is a clumping factor dependent on the vegetation class (Anderson et al., 2005). TS and TC are used to compute the surface energy budget for the soil and canopy components separately:

RNS = HS + LES + G

(8)

RNC = HC + LEC

(9)

where RNS is net radiation at the soil surface, RNC is net radiation divergence in the canopy layer, LES is the soil evaporation rate, LEC is the canopy transpiration rate, HS and HC are the soil and canopy sensible heat fluxes, and G is the soil heat flux. Considering a network of resistances in series (Norman et al., 1995; Kustas and Norman, 1999), HS and HC can be expressed as:

HS = ρCP

TS − TAC RS

(10)

HC = ρCP

TC − TAC RX

(11)

3.3. Plot scale TSEB input data The TSEB model was applied at the plot scale, with local meteorological input, at the two eddy covariance stations established at the irrigated meadow and at the rainfed pasture in Val Mazia/Matschertal. As input data for the plot scale simulations, we used a combination of ground and satellite-based observations. Specifically, the radiometric temperature, TRAD, was calculated from the upwelling longwave component of a 4-component net radiometer, shortwave radiation was measured with a pyranometer, net radiation with the net radiometer, and meteorological parameters were measured by the local meteorological stations. Daily LAI was retrieved by a spline interpolation of the 4 days-1 km MODIS standard product (Cammalleri et al., 2014a). Vegetation height, hc, was calculated according to Eq. 10 of Anderson et al. (2007):

with the total H corresponding to:

H = ρCP

TAC − TA RA

(12)

where ρ is air density, CP is the heat capacity of air, TAC is the temperature of the air within the canopy, RX is the total boundary layer resistance of the canopy, RS is the soil resistance to the transport of sensible heat, and RA is the aerodynamic resistance (Kustas and Norman, 1999; Kustas and Norman, 2000a,b). The canopy component of the latent heat flux, LEC, is calculated by the Priestley-Taylor(PT) equation (Priestley and Taylor, 1972) applied to the divergence of net radiation within the canopy:

LEC = αPTC fG

Δ RNC Δ+γ

hc = h min + f (0)[h max − h min]

where f(0) is the vegetation cover fraction for sensor at nadir, which has been derived from field observations, and hmin and hmax are the seasonal nominal minimum and maximum canopy heights for the land-cover class grassland. The fraction of vegetation in canopy that is green was derived from MODIS normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) at 250 m resolution, after Fisher et al. (2008):

(13)

where αPTC is a coefficient related to the Priestley-Taylor coefficient but derived only for the canopy, Δ is the slope of the saturation vapor pressure versus temperature curve, fG is the fraction of green vegetation, and γ is the psychrometric constant. A first guess value of LEC is calculated with an initial value of αPTC, set to 1.26 (Agam et al., 2010), which is iteratively reduced under stress conditions. The latent heat flux from the soil surface is calculated as the residual of the energy balance equation:

LES = RNS − G − HS .

fg = 1.2

EVI . NDVI

(17)

3.4. ALEXI/DisALEXI input data The ALEXI simulations used in this study exploited as input LST, leaf area index (LAI), downwelling shortwave and longwave radiation generated by the LSA-SAF (Land Surface Analysis Satellite Application Facility, http://landsaf.meteo.pt/) from Meteosat Second Generation (MSG) satellite data. The MSG nominal spatial resolution over Italy is of approximately 3.5 km. The meteorological variables (air temperature, wind speed, vapor pressure) required by ALEXI were simulated by the

(14)

G is estimated as a fraction of RNS:

G = cG RNS

(16)

(15) 330

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

Table 1 R2, MBD [W m−2], and RMSE [W m−2] between daytime (6–18 UTC+1) hourly averages of observed and TSEB modeled latent (LE) and sensible (H) heat at the meadow site for the years 2012, 2013 and 2014, including data from March to November. h, le, rn and g are the averages of the observed fluxes. Year

2012 2012 2012 2013 2013 2013 2014 2014 2014

Closure

Res BR NC Res BR NC Res BR NC

H

LE

RN

G

h

R2

MBD

RMSE

le

R2

MBD

RMSE

rn

R2

MBD

RMSE

g

R2

MBD

RMSE

53 123 53 57 130 57 64 157 64

0.85 0.63 0.85 0.69 0.57 0.69 0.66 0.35 0.66

0.5 −51 0.5 −7 −64 −7 −23 −99 −23

25 84 25 37 100 37 45 139 45

222 202 127 224 202 128 226 198 115

0.96 0.96 0.81 0.95 0.94 0.77 0.94 0.94 0.75

−4 14 46 0.1 19 52 2 29 68

28 34 80 31 40 88 32 45 101

270 – – 273 – – 226 – –

0.98 – – 0.99 – – 0.97 – –

15 – – 13 – – 23 – –

32 – – 23 – – 43 – –

16 – – 17 – – 9 – –

0.58 – – 0.53 – – 0.14 – –

29 – – 24 – – 41 – –

29 – – 38 – – 61 – –

model ALEXI/DisALEXI in the study area. We evaluated hourly averages of the TSEB modeled sensible and latent heat fluxes exploiting EC measurements. Modeled net radiation was compared with net radiometer measurements, and modeled soil heat flux was compared with soil heat plate measurements corrected for the heat storage between the soil surface and the instrument. We confronted TSEB simulations and EC fluxes at the meadow station, from 2012 to 2014 (Table 1, Fig. 3 (a)–(i)), and at the pasture site for years 2014 and 2015 (Table 2, Fig. 4 (a)–(f)). Results strongly depend on the method used for forcing the closure of the energy balance. In the meadow, RMSE for H and LE is between 25 and 45 W m−2 with Res, between 34 and 139 W m−2 with BR, and between 25 and 101 W m−2 with NC. For RN, RMSE is between 23 and 43 W m−2, and for G it is between 29 and 61 Wm−2. In the pasture, RMSE for H and LE ranges between 42 and 73 W m−2 with Res, 31 and 114 W m−2 with BR, and 57 and 105 W m−2 with NC. RMSE for RN is of 16–19 W m−2 and for G it is 61–66 W m−2. In conclusion, at the meadow site, fluxes corrected by the LE-residual method agree better with TSEB simulations than uncorrected fluxes or fluxes corrected by the BR method. This outcome is consistent with the assumption of validity for our study site of the results of Wohlfahrt et al. (2010), who compared EC measurements with independent measurements of ET, and showed that the Res method performs better than the BR method in temperate mountain grasslands. In contrast, at the pasture site, fluxes corrected by the BR method agree better with TSEB simulations than fluxes corrected by the Res method.

Weather Research and Forecast model (WRF), with an output frequency of 30 min and a spatial resolution of 3 km. For the DisALEXI simulations, we used the following remote sensing input: MODIS LAI (MCD15A3), NDVI (MOD13A2), albedo (MCD43B3), geolocation (MOD03) and LST (MOD11 L2), Landsat ETM+ thermal and optical imagery (the only one available in 2012). High resolution land-cover maps of the Provinces of Bolzano/Bozen and Trento were resampled at MODIS and Landsat resolution. Before running DisALEXI, MODIS and Landsat data were preprocessed. First Landsat thermal imagery was atmospherically corrected by the single channel method described in Jiménez-Muñoz and Sobrino (2003) and Sobrino et al. (2004), as revised in Jiménez-Muñoz et al. (2009), using MODIS water vapor product, MOD05. Then Landsat LAI was generated with a regression tree approach trained with samples from the MODIS LAI products and Landsat reflectances (Gao et al., 2012). Regarding MODIS, smoothed and gapfilled daily LAI maps were generated by the TIMESAT algorithm (Jönsson and Eklundh, 2004) from the 4-day MODIS LAI product, MCD15A3. 4. Results and discussion 4.1. TSEB simulations at the plot scale The TSEB was first applied to local meteorological data collected in proximity of the EC towers. These simulations were considered as benchmark for evaluating the performances of the remote-sensing

Fig. 2. Cumulative flux footprint of the meadow site calculated based on Kljun et al. (2015), with a Bing map on the background. The isolines of 10 to 90% flux footprint contributions are indicated by red lines, and the flux tower is indicated by a black dot. The picture was taken during a moment when the meadow was freshly cut. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

331

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

(a) Residual closure - year 2012

(b) Bowen ratio preservation

(c) No closure

(d) Residual closure - year 2013

(e) Bowen ratio preservation

(f) No closure

(g) Residual closure - year 2014

(h) Bowen ratio preservation

(i) No closure

Fig. 3. Scatterplots of modeled vs observed fluxes at the eddy covariance station in the meadow in Val Mazia/Matschertal for the years 2012–2014. Measured fluxes have been corrected as indicated below each figure.

It was not possible to evaluate Landsat based fluxes because at the time of Landsat overpass (around 10 UTC) the pixels corresponding to the EC grassland measurement sites in Val Mazia/Matschertal were always covered by clouds. This problem is probably due to the influence of complex topography on the formation of convective clouds, which causes overcast conditions at the time of Landsat overpass (Fig. 5). The occurrence of clouds in combination with the stripes present in ETM+

4.2. DisALEXI simulations at the landscape scale 4.2.1. Modeled fluxes evaluation at the overpass time DisALEXI simulations driven by MODIS and Landsat data were performed over South Tyrol for the years 2012 and 2013. Modeled fluxes at the time of satellite overpass were compared with EC measurements.

Table 2 R2, MBD [W m−2], and RMSE [W m−2] between daytime (6–18 UTC+1) hourly averages of observed and TSEB modeled latent (LE) and sensible (H) heat at the pasture site for 2014 (from June to December), and for 2015 (from March to November). h, le, rn and g are the averages of the observed fluxes. Year

2014 2014 2014 2015 2015 2015

Closure

Res BR NC Res BR NC

H

LE

h

R

2

MBD

RMSE

84 102 84 132 139 132

0.85 0.66 0.85 0.83 0.48 0.83

19 37 19 −40 62 −40

57 83 57 73 114 73

RN

G

le

2

R

MBD

RMSE

rn

R

MBD

RMSE

g

R2

MBD

RMSE

168 156 133 182 172 154

0.86 0.93 0.66 0.84 0.89 0.42

−8 −8 −26 −50 −46 −61

42 31 71 66 57 105

199 – – 236 – –

1 – – 0.99 – –

-3 – – 5 – –

16 – – 19 – –

32 – – 40 – –

0.1 – – 0.1 – –

−26 – – −19 – –

61 – – 66 – –

332

2

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

(a) Residual closure - year 2014

(d) Residual closure - year 2015

(b) Bowen ratio preservation

(c) No closure

(e) Bowen ratio preservation

(f) No closure

Fig. 4. Scatterplots of modeled vs observed fluxes at the eddy covariance station in the pasture in Val Mazia/Matschertal for the years 2014–2015. Measured fluxes have been corrected as indicated below each figure.

imagery, generate gaps which are too wide to be filled. Direct comparisons between tower observations and DisALEXI flux retrievals using MODIS data at 1 km resolution are difficult due to mismatch between the model pixel and tower footprint scales. An estimate of the footprint extension for the meadow EC site in Val Mazia/ Matschertal is shown in Fig. 2, with 90% of fluxes originated from an area homogeneously covered by grassland, extending 50 m north and 100 m south, and 50 m east and west of the flux tower. In contrast, the modeled response incorporates different land-cover types falling in the same 1-km pixel (forest, shrubs, as in Fig. 1). In addition, due to the fragmented ownership of the fields, the small patches which are included in a MODIS pixel are not irrigated simultaneously, and have different soil moisture conditions. While direct assessment of model accuracy is not possible given this level of scale mismatch, we can still use the tower fluxes to evaluate the relative partitioning of energy flux components as modeled by DisALEXI. Scatterplots comparing observed and modeled energy fluxes at the time of MODIS overpass in 2012 and 2013 are represented in Fig. 6 (a), (b) and (c). In these comparisons, EC measurements were corrected using the Res method (Fig. 6 (a)), the BR method (Fig. 6 (b)), and without correction (Fig. 6 (c)). Statistical metrics of model-measurement agreement for the years 2012 and 2013, including MBD and RMSE, are summarized in Table 3. Model discrepancies for RN and G are comparable with those from the plot scale simulations (Section 4.1), with an RMSE of 54–76 W m−2 for RN and 17–32 W m−2 for G. In contrast, discrepancies for H and LE are much higher, and strongly influenced by the method used for forcing the energy balance closure. RMSEs of 109–154 W m−2 and 100–172 W m−2 are obtained for H and LE, respectively, with the lower values associated with BR closure and

higher values with Res. Tower LE is underestimated in all cases at the model pixel scale, with MBD between −76 W m−2 and −143 W m−2 for closure enforced fluxes. This consistent bias may be related to the fact that the area contributing to the tower fetch is irrigated, while the 1 km model pixel also includes unirrigated land. These discrepancies are higher than obtained in previous TSEB simulations based on MODIS data performed on grassland ecosystems in flat regions. For example, Guzinski et al. (2013) report RMSE between 76 and 89 W m−2 for H, and between 120 and 128 W m−2 for H, for semi-arid and temperate grassland in Denmark, while for RN and G results are similar to ours. In addition, Gan and Gao (2015) obtained a RMSE of 63 W m−2 for H and 94 W m−2 for LE in a grassland in the SongNen plain of China. However, these sites were more homogeneous at the 1 km scale than is the Val Mazia site. In addition to the scale mismatch between the model pixel and tower footprint, the uncertainties of the input data due to the retrieval process cause additional discrepancy from EC measurements. In the following section, we perform a quantitative analysis on the influence of the uncertainty of the input data on model performance. 4.2.2. TSEB sensitivity to LAI, LST and shortwave radiation In this study, we quantified the model error attributable to the uncertainty of the remote sensing input data by performing a sensitivity study on the TSEB model at the plot scale. We specifically focused on LST, LAI and solar radiation (shortwave radiation, SWR). The reason behind the choice of these variables is that they are particularly influenced by the heterogeneity of the study area in terms of land-cover and topography (altitude, slope, exposition). In fact, the spatial variability of LST and LAI is dominated by land-cover, while solar radiation is

333

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

Fig. 5. LST over the Province of Bolzano/Bozen derived from Landsat thermal data in the year 2012. The figure shows all the available Landsat ETM+ images which were less affected by cloud coverage. An asterisk is placed at the measurement site of Val Mazia/Matschertal.

found a RMSD of 130 W m −2 between observed incoming solar radiation and the LSA-SAF product at the time of MODIS overpass (Fig. 7). The sensitivity was calculated in W m −2 as the absolute difference between the energy flux value modeled by decreasing or increasing the variable X of δX, and the original value. For example, the sensitivity of H to an increase of LST was calculated as:

strongly influenced by topography. Since the resolution of the satellite products cannot capture small scale heterogeneities, we can expect significant differences between the pixel values and the local true values of LST, LAI and SWR. For running DisALEXI, we exploited standard MODIS LST and LAI, with a resolution of 1 km. The expected accuracy of LST is within 1 K for areas with known emissivities in the range −10 °C to 50 °C (Wan et al., 2004), while the accuracy of LAI (MCD15 product, Terra+ Aqua) for herbaceous biomes is 1.16 (Fang et al., 2012). In contrast, SWR is the operational LSA product based on SEVIRI data with a nominal resolution of 3.5 km. Validation studies show a RMSE of 90–120 W m −2 (LSA-SAF, 2011). Consequently, in the sensitivity analysis we varied LST by ± 1 K, LAI by ± 1.16, and SWR by ± 100 W m −2. However, in our study area we expect even higher differences between ground measured and satellite retrieved solar radiation, due to the small scale variability of topography and vegetation. In fact, for the year 2012 we

δH = |H (LST + δLST ) − H (LST )|

(18)

where δH is the absolute difference between H obtained increasing LST of δLST, and the original H value. For better approximating the effect of errors in LST, LAI and SWR on fluxes modeled by DisALEXI, we considered only the time interval 9 a.m.–12 a.m. UTC, which includes all MODIS daytime overpasses, because the model sensitivity exhibits a diurnal cycle, as a consequence of the diurnal cycle of the available energy.

334

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

(a) Resisual closure

(b) Bowen ratio preservation

(c) No closure Fig. 6. Energy fluxes observed and modeled by DisALEXI driven by MODIS data in 2012 and 2013 for the meadow site in Val Mazia/Matschertal at the time of MODIS overpass. Eddy covariance fluxes have been corrected in Fig. 6(a) by attributing all the residual to LE, in Fig. 6(b) by preserving the Bowen ratio, and have not been corrected in Fig. 6(c).

SWR. The influence of LST variation of ± 1 K was generally below 20 W m −2. The effect of LAI variation was higher for H than for LE, with values between 10 and 50 W m −2 for H, decreasing with increasing LAI, and below 20 W m −2 for LE. In summary, the most critical variable is solar radiation, whose uncertainty can cause a variation of LE larger than 100 W m −2, thus

The results of the sensitivity analysis performed on the TSEB model are represented in Figs. 8 and 9 as boxplots of δH and δLE for an increase and a decrease of LST, LAI and SWR. The alteration of SWR produced a strong variation in LE, with median values between 60 and 70 W m −2. We found a big spread of δLE for an increase of SWR, reaching values above 100 W m −2, and a small spread for a decrease of

Table 3 MBD [W m−2], and RMSE [W m−2] for energy fluxes modeled by DisALEXI at the meadow station in Val Mazia/Matschertal at the time of MODIS overpass for 2012 and 2013, from April to October. In brackets RMSE is expressed as percentage of the mean observed value. Results are given correcting EC observations of LE and H with the residual method (Res) or preserving the Bowen ratio (BR). Year

2012 2012 2012 2013 2013 2013

Closure

Res BR NC Res BR NC

H

LE

RN

G

MBD

RMSE

MBD

RMSE

MBD

RMSE

MBD

RMSE

142 80 139 126 66 126

148 109 143 154 118 154

−108 −76 28 −143 −76 27

179 142 131 166 117 100

−33 – – −60 – –

54 (9%) – – 76 (12%) – –

11 – – 11 – –

17 (54%) – – 32 (87%) – –

(168%) (83%) (180%) (121%) (64%) (121%)

335

(47%) (41%) (54%) (40%) (32%) (40%)

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

56 PM) reference ET for a hypothetical grass, as described by Allen et al. (1998). This indicator is called ESI (Evaporative Stress Index). The ESI has been compared with other standard drought indexes, showing a good correspondence with short-term drought metrics like the National Climate Data Center (NCDC) Palmer Z index (Anderson et al., 2011). The advantage of the ESI is that it does not require precipitation measurements and it can be produced at higher resolution than standard drought indexes. Given the limited temporal extent of the analyses presented here, in this study we evaluate fPET map timeseries rather than anomalies in fPET. Daily ET was calculated by upscaling the instantaneous LE estimated by DisALEXI, relying on the assumption that the evaporative fraction is constant during the day (Brutsaert and Sugita, 1992):

ETd =

LE (t ) 1.1 × (RNd − Gd ) × [mm day −1] RN (t ) − G (t ) 2.45 −2

where all the energy fluxes are expressed in MJ m day , (RNd − Gd) is the available energy from sunrise to sunset, the coefficient 2.45 MJ kg−1 is the latent heat of vaporization of water at 20 °C, and the constant 1.1 accounts for the daytime underestimation of ET and for the omission of the nighttime ET. After upscaling to daily ET, a gapfilling method was applied based on conservation and smoothing the ratio of actual-to-reference ET (Anderson et al., 2007). This procedure mitigated the impact of cloud cover which strongly affects data availability in the study area. A comparison between DisALEXI and measured daily ET would have been necessary to assess the accuracy of daily ET, but no days were found without gaps in the LE measured by the eddy covariance system. Consequently, for a rough indication of the representativeness of daily ET, we compared DisALEXI ET with TSEB ET. We calculated TSEB ET by cumulating the complete time series of hourly LE. Then we quantified the error due to the upscaling of daily fluxes from instantaneous model output. We first compared the TSEB cumulated ET with the TSEB ET upscaled by Eq. (20) (Fig. 10), then we compared DisALEXI and TSEB upscaled ET (Fig. 11), and finally DisALEXI ET with TSEB cumulated ET (Fig. 12). We calculated the ratio of RMSE between TSEB cumulated ET and TSEB scaled ET, and RMSE between MODIS scaled

Fig. 7. Comparison between observed incoming solar radiation and LSA-SAF solar radiation at the time of MODIS overpass in 2012 at the measurement station in the meadow, in Val Mazia/Matschertal.

strongly influencing the accuracy of the model output. 4.2.3. Daily ET and evaporative stress Anderson et al. (2008) and Anderson et al. (2016), among other works, show the potential of the drought indicator defined as the standardized anomaly in the ratio between actual and potential evapotranspiration (PET), fPET:

fPET =

ET PET

(20) −1

(19)

where ET is the daily actual ET, and PET is the Penman-Monteith (FAO-

Fig. 8. Boxplots showing the sensitivity of H calculated by the TSEB model to an increase and a decrease of LST, LAI and SWR, between 9 a.m. and 12 a.m. UTC for year 2012. One outlier was excluded from the plot for graphical reasons. Specifically, the plot does not include the sensitivity of H to an increase of SWR obtained on the 18th of May at 11:00, equal to 270 W m−2).

Fig. 9. Boxplots showing the sensitivity of LE simulated by the TSEB model to an increase and a decrease of LST, LAI and SWR, between 9 a.m. and 12 a.m. UTC for year 2012. Two outliers were excluded from the plot for graphical reasons. Specifically, the plot does not include the sensitivity of LE to an increase of SWR obtained on the 18th of May at 11:00 and on the 3rd of August at 13:00, equal respectively to to 373 and 478 W m−2).

336

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

Fig. 10. Daily ET in the year 2012 calculated by cumulating hourly TSEB LE and by upscaling instantaneous TSEB LE assuming the evaporative fraction constant during the day (Eq. (20)).

Fig. 11. Comparison between DisALEXI and TSEB ET in the year 2012 upscaled from instantaneous midday LE assuming the evaporative fraction constant during the day (Eq. (20)).

Fig. 12. Daily ET [mm day−1] in the year 2012 upscaled by DisALEXI assuming the evaporative fraction constant during the day (Eq. (20)), and cumulated by TSEB from local hourly data.

337

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

Fig. 13. fPET in the Provinces of Bolzano/Bozen (North) and Trento (South) derived applying DisALEXI to MODIS data.

analyze, also because we do not know the amount of irrigation. At Rovereto, the fPET timeseries follows the trend of CSI very well, both showing maximum of stress (minimum of fPET, maximum of CSI) in August during a period of low rainfall. At Passo Rolle there is abundant rain and low potential evapotranspiration due to the high elevation of the site (2000 m a.s.l.), which is in agreement with the absence of evaporative stress and with low values of CSI. At Silandro, a waterlimited site located in the North-Eastern side of the Province of Bolzano/Bozen, both CSI and fPET are high during summer. This evidence could be attributed to irrigation. In fact, CSI is based on the natural precipitation measured in this location, while our fPET estimation is based on remote sensing observations, which are able to capture the effect of irrigation on energy partitioning. For a further evaluation of the fPET signals, we consulted maps of Combined Drought Index (CDI) generated by the European Drought Observatory of the Joint Research Centre EDO JRC, http://edo.jrc.ec. europa.eu), which combines the Standardized Precipitation Index (SPI), the soil moisture anomalies, and the Fraction of Absorbed Photosynthetically Active Radiation (fAPAR) anomalies. Despite the low resolution (0.0416 deg), the CDI maps confirmed the presence of vegetation stress in September in the Adige Valley(where station 4 of

ET and TSEB cumulated ET. From this ratio, we estimated that 22% of the RMSE between MODIS ET and TSEB ET can be roughly attributed to the upscaling method. Once compared ET with point simulations, daily maps of fPET were developed over the study area using daily upscaled and gapfilled MODIS DisALEXI ET and gridded PET computed using the meteorological input dataset described in Section 3.4. fPET maps for June to September of 2012 over the Provinces of Trento and Bolzano/Bozen are shown in Fig. 13, focusing on maps that were minimally impacted by residual cloud cover. These maps show the increasing of evaporative stress from June to September, especially in Val Venosta (North-West area in Fig. 13, including station 3 of Fig. 1) in the province of Bolzano/ Bozen, and Val d’Adige (South-central area in Fig. 13, including Rovereto, the station 4 of Fig. 1) in the Province of Trento. As a comparative metric, we computed the crop stress index (CSI) at six stations located in grassland in the Provinces of Trento and Bolzano/ Bozen (Fig. 14), whose position is shown in Fig. 1. The CSI was calculated according to Gage et al. (2015), as the ratio between accumulated temperature and accumulated precipitation, summed over a month long period. All the stations are in irrigated grasslands, a part from Rovereto and Passo Rolle, the latter being the most interesting to 338

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

Fig. 14. fPET and crop stress index (CSI) for the year 2012 at six measurement stations located in grasslands in the Provinces of Trento and Bolzano/Bozen (precipitation and air temperature data from http://www.meteotrentino.it/ and http://www.provincia.bz.it/meteo/home.asp). The position of the stations is showed in Fig. 1

compensated by irrigation. One reason for the agreement between DisALEXI and measurements is the favorable situation that, despite the low resolution of MODIS, irrigated meadows are the predominant landcover type in the pixel including the EC station in the meadow, and rainfed pastures cover most of the pixel including the EC station in the pasture.

Fig. 14 is located), and also the early stress in June at East of the Garda Lake, which is also evident in Fig. 13 at South-East of the Province of Trento. As an additional test of MODIS-based fPET utility in this study area, we investigated whether this product could capture the contrast of water availability conditions between irrigated meadows and pastures that are not irrigated in Val Mazia /Matschertal. To this end, we compared the temporal patterns of daily fPET from MODIS DisALEXI with precipitation and CSI computed from local meteorological data collected at the two flux measurement sites (Fig. 15). Irrigation was not measured at these sites, but since the two stations are very close, we assumed that the difference in the measured precipitation was equivalent to irrigation. Despite the coarse resolution of the MODISdriven DisALEXI retrievals, we found that the model is capable of identifying the contrasting conditions of these two land-cover types. In fact, starting from August the pasture site experienced a decrease in fPET, which corresponded to very low precipitation levels from the end of June, and to high values of CSI from the middle of July to September (Fig. 15 shows the patterns of ET, fPET, the 30-days cumulative precipitation, and CSI). The meadow shows very low levels of stress compared to the pasture, because the lack of precipitation was

5. Summary and conclusions In this study, we applied the model ALEXI/DisALEXI for the first time in a mountainous region, specifically in Trentino-Alto Adige, in Italy. We especially focused on Val Mazia/Matschertal, in the Northwest of the Province of Bolzano/Bozen, which is one of the driest areas in the Alps (Isotta et al., 2014). Two eddy covariance towers in managed and unmanaged grasslands were used for evaluating the model. Results demonstrate that the model is affected by big limitations, especially due to the lack of high resolution thermal data in the area of interest. In fact, the main input of ALEXI/DisALEXI is radiometric temperature, which we derived from MODIS and Landsat polar orbiting satellites. Unfortunately, only MODIS-based retrievals of energy fluxes could be analyzed and evaluated in the area of interest, because during

339

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

remote sensing input data of solar radiation and leaf area index. Since one known issue in complex terrain is the influence of topography on incoming solar radiation, and consequently on net radiation, due to small-scale changes in slope, aspect and elevation, we hypothesize that DisALEXI would benefit of higher resolution solar radiation data, like the one produced by the HelioMont algorithm (Stöckli, 2013), whose hight accuracy in the Alps has been assessed by Castelli et al. (2014). Another issue in the area of interest is the small scale variability of landcover, and consequently of leaf area index (LAI) and vegetation height. As evident from the sensitivity analysis, LAI is also very influential on DisALEXI output, and its accuracy could be improved. For example, concerning grasslands, Pasolli et al. (2015a) show that the standard MODIS LAI product is affected by several issues in mountainous regions, and propose an alternative algorithm for retrieving LAI from MODIS data, which is specifically tuned to mountain grasslands. While instantaneous fluxes were directly compared with EC data, the upscaled daily ET was confronted with the cumulative ET calculated from hourly TSEB simulations, because EC data were affected by gaps. An RMSE of 1.89 mm day−1 was found between DisALEXI upscaled and TSEB cumulated ET in Val Mazia/Matschertal in the year 2012, which is partially attributable to the upscaling method (Cammalleri et al., 2014b). The upscaled ET was used for mapping a thermal stress index, using the ratio of actual-to-potential ET (fPET). We compared the fPET with cumulative precipitation for different grassland stations in the area of interest, and observed a good correspondence. The fPET was also in agreement with the drought index CDI. In addition, despite its moderate resolution, we found that MODIS-based fPET was able to catch clearly the differences between managed and unmanaged grasslands in the study area. This suggests that the fPET can also be exploited in mountainous regions for monitoring water stress and droughts at the regional scale. The possibility of retrieving this index is quite relevant, because it reveals water stress without distributed precipitation measurements and without information of the soil hydraulic properties, which are hard to obtain in complex terrain, with high spatial variability of soil type. As stated in the Introduction, this work was carried out simultaneously with the estimation of soil moisture from microwave remote sensing data (Pasolli et al., 2015b). We foresee to extend the analysis by combining the strengths of the two approaches for monitoring water stress conditions. From one side, in fact, methods based on thermal data can provide data over dense vegetation, and from the other side microwave methods are also available under cloudy sky conditions. One possibility is to derive a thermal-based soil moisture, and assimilating it into a water balance model (Hain et al., 2012), together with microwave-based soil moisture. Another one is calculating a water stress index by integrating a representation both of the root-zone conditions, diagnosed by the thermal based ET, and of the surface conditions, investigated by microwave remote sensing. In conclusion, DisALEXI has a large potential for application in mountainous areas, such as the Alps, and there is room for improving the accuracy of the remote sensing ancillary data, but a major limitation comes from the scarcity of high resolution thermal data. For thermal remote sensing approaches to have practical utility for modeling daily turbulent fluxes in mountainous regions with persistent cloud cover and complex terrain, the temporal acquisition frequency of high resolution thermal remote sensing must be significantly increased. Planned missions like the NASA-JPL HyspIRI (Hyperspectral Infrared Imager) (Hochberg et al., 2015), and the JPL ECOSTRESS (The ECOsystem Spaceborne Thermal Radiometer Experiment on Space Station) (Guillevic et al., 2014) offer good examples of thermal imagery which could strongly support ET mapping in areas characterized by complex terrain and heterogeneous landscape, by combining high spatial resolution and short revisit time. Specifically, the HyspIRI mission, whose launch is foreseen after 2022, includes a multispectral imager with 8 spectral bands between 4 and 12 μ m, a spatial resolution of 60 m at

Fig. 15. MODIS-based ET and fPET, 30 days cumulative rain, and crop stress index (CSI) in the year 2012 at the meadow and pasture measurement stations in Val Mazia/ Matschertal.

the period of investigation Landsat pixels were always contaminated by clouds. While we expect similar problems for all models that exploit thermal data in complex terrain, the TSEB approach has the advantage of being process-oriented and capable of separating soil and canopy contributions to ET. Before applying the remote sensing model, we performed plot scale TSEB model simulations driven by local data, and we compared the results of these simulations with eddy covariance data. We found that the TSEB model reproduces LE with a RMSE between 28 and 104 W m −2 , depending on the method that is used to correct EC fluxes for forcing the energy budget closure. In particular, in the meadow we obtained the best agreement between EC and TSEB fluxes by attributing all the residual of the energy budget to LE, whereas in the pasture we obtained the best agreement preserving the Bowen ratio. Since the lack of closure in the energy budget is an important issue in our mountain stations, additional effort is needed for understanding and quantifying the causes of such an unclosure, so that EC fluxes can be a more robust reference to validate energy balance models. The plot scale simulations were considered as a benchmark for the modeling based on remote sensing. MODIS-based DisALEXI 1 km simulations were performed covering the entire region Trentino-Alto Adige and evaluated against EC data in terms of instantaneous fluxes estimated at the time of satellite overpass. Results evidence large discrepancies between modeled and observed LE and H, exceeding 100 W m −2. Much of this error may be due to significant scale mismatch between the tower flux footprint (on the order of 100 m) and the 1 km model pixel. A sensitivity analysis performed on TSEB demonstrated that additional sources of error may be related to inaccuracies in 340

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

LAI LEC LES LST LTER MBD METRIC MODIS MSG MW NDVI RMSE R RA RN RNS RNC RS RX SEBAL SPI SWR TA TAC TC θ TIR TRAD TS TSEB VIIRS WMO WRF Ω

nadir, a swath of 600 km, and a revisit of 5 days, while ECOSTRESS, whose launch is foreseen in 2017, will use the PHyTIR Thermal Infrared Radiometer, with 5 spectral bands between 8 and 12 μ m, a spatial resolution of 68 m × 38 m (crosstrack × downtrack), a swath of 400 km, and a revisit of 4 days. It is also worth noting that when two fully functional Landsat satellites will be available, upon launch of Landsat 9 in 2020, the opportunities for clear-sky retrievals at 100 m will be improved, and that there is discussion going on about adding Landsat-resolution thermal capabilities to the Sentinel series (Sentinel 8). Moreover, some ongoing modifications in the modeling system may significantly improve DisALEXI performance over mountainous landscape. These modifications include: i) the adaptation of the package to MODIS Collection 6 products, which will allow running the MODIS disaggregation at 500 m rather than 1 km, and ii) a new ALEXI approach using day-night temperature differences obtained with the new VIIRS instrument which has a 375 m thermal infrared channel and provides MODIS-like near daily coverage. In addition to the need for high resolution thermal remote sensing data, suitably gridded ground observations would be necessary for a more robust evaluation of ALEXI/DisALEXI, and for improving energy flux modeling based on remote sensing data in complex terrain. There are several FLUXNET stations in the Alps which are measuring the turbulent fluxes and all the meteorological variables, and it has been proposed to add LAI measurements. There are also a few FLUXNET sites in the mountain areas of the United States and in Tibet; however, these sites are sparsely distributed. A further effort in the direction of a properly equipped mountain area will be the expansion of the instrumental network of the Long Term (Socio) Ecological Research area Mazia/Matsch in the Italian Alps. Acknowledgments The work presented in this paper was supported by the projects HiResAlp and MONALISA funded by the Province of Bolzano/Bozen within the Legge 14 research program. The authors would like to thank the MODIS and Landsat teams for providing optical and thermal imagery, and the meteorological services of the Provinces of Trento and Bolzano/Bozen for providing access to historical weather data.

Leaf Area Index canopy latent heat flux coil latent heat flux land surface temperature Long Term Ecological Research mean bias deviation Mapping Evapotranspiration with Internalized Calibration MODerate resolution Imaging Spectroradiometer Meteosat Second Generation microwave Normalized Difference Vegetation Index Root Mean Square Error linear regression coefficient aerodynamic resistance net radiation net radiation at the soil surface net radiation divergence in the canopy layer soil resistance to the transport of sensible heat boundary layer resistance of the canopy Surface Energy Balance Algorithm for Land Standardized Precipitation Index shortwave radiation temperature of the air above the canopy temperature of the air within the canopy canopy temperature thermal sensor view angle thermal infrared radiometric temperature soil temperature two source energy balance Visible Infrared Imaging Radiometer Suite World Meteorological Organization weather research and forecast clumping factor

References Agam, N., Kustas, W.P., Anderson, M.C., Norman, J.M., Colaizzi, P.D., Howell, T.A., Prueger, J.H., Meyers, T.P., Wilson, T.B., 2010. Application of the Priestley-Taylor approach in a two-source surface energy balance model. J. Hydrometeorol. 11, 185–198. Allen, R.G., Pereira, L.S., Raes, D., Smith, M., et al., 1998. Crop evapotranspirationguidelines for computing crop water requirements — FAO irrigation and drainage paper 56. FAO, Rome 300 D05109. Allen, R.G., Tasumi, M., Trezza, R., 2007. Satellite-based energy balance for mapping evapotranspiration with internalized calibration (METRIC)-model. J. Irrig. Drain. Eng. 133, 380–394. Anderson, M., Norman, J., Diak, G., Kustas, W., Mecikalski, J., 1997. A two-source timeintegrated model for estimating surface fluxes using thermal infrared remote sensing. Remote Sens. Environ. 60, 195–216. Anderson, M., Norman, J., Kustas, W., Houborg, R., Starks, P., Agam, N., 2008. A thermalbased remote sensing technique for routine mapping of land-surface carbon, water and energy fluxes from field to regional scales. Remote Sens. Environ. 112, 4227–4241. Anderson, M.C., Hain, C., Wardlow, B., Pimstein, A., Mecikalski, J.R., Kustas, W.P., 2011. Evaluation of drought indices based on thermal remote sensing of evapotranspiration over the continental United States. J. Climate 24, 2025–2044. Anderson, M.C., Kustas, W.P., Alfieri, J.G., Gao, F., Hain, C., Prueger, J.H., Evett, S., Colaizzi, P., Howell, T., Chávez, J.L., 2012. Mapping daily evapotranspiration at landsat spatial scales during the BEAREX’08 field campaign. Adv. Water Resour. 50, 162–177. Anderson, M.C., Norman, J., Kustas, W.P., Li, F., Prueger, J.H., Mecikalski, J.R., 2005. Effects of vegetation clumping on two-source model estimates of surface energy fluxes from an agricultural landscape during smacex. J. Hydrometeorol. 6, 892–909. Anderson, M.C., Norman, J.M., Mecikalski, J.R., Otkin, J.A., Kustas, W.P., 2007. A climatological study of evapotranspiration and moisture stress across the continental united states based on thermal remote sensing: 1. Model formulation. J. Geophys. Res.: Atmos. (1984-2012) 112. Anderson, M.C., Zolin, C.A., Sentelhas, P.C., Hain, C.R., Semmens, K., Yilmaz, M.T., Gao, F., Otkin, J.A., Tetrault, R., 2016. The evaporative stress index as an indicator of agricultural drought in Brazil: an assessment based on crop yield impacts. Remote Sens. Environ. 174, 82–99. http://www.sciencedirect.com/science/article/pii/ S0034425715302212http://dx.doi.org/10.1016/j.rse.2015.11.034. Aubinet, M., Grelle, A., Ibrom, A., Rannik, Ü., Moncrieff, J., Foken, T., Kowalski, A.,

List of acronyms and symbols αPT ABL ALEXI BR CDI Δ

Priestley-Taylor coefficient Atmospheric Boundary Layer Atmosphere Land EXchange Inverse Bowen ratio Combined Drought Index slope of the saturation vapor pressure versus temperature curve DEM Digital Elevation Model DisALEXI Disaggregate ALEXI EC eddy covariance EDOJRC European Drought Observatory of the Joint Research Centre EEA European Environment Agency ET evapotranspiration ETM+ Enhanced Thematic Mapper Plus EVI Enhanced Vegetation Index f vegetation cover fraction fAPAR fraction of Absorbed Photosynthetically Active Radiation fG fraction of green vegetation γ psychrometric constant G soil heat flux hc canopy height HC canopy sensible heat flux HS soil sensible heat flux 341

Remote Sensing of Environment 209 (2018) 327–342

M. Castelli et al.

(HyspIRI). Isotta, F.A., Frei, C., Weilguni, V., Perčec Tadić, M., Lassègues, P., Rudolf, B., Pavan, V., Cacciamani, C., Antolini, G., Ratto, S.M., et al., 2014. The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine rain-gauge data. Int. J. Climatol. 34, 1657–1675. Jiménez-Muñoz, J.C., Cristóbal, J., Sobrino, J., Soria, G., Ninyerola, M., Pons, X., et al., 2009. Revision of the single-channel algorithm for land surface temperature retrieval from landsat thermal-infrared data. IEEE Trans. Geosci. Remote Sens. 47, 339–349. Jiménez-Muñoz, J.C., Sobrino, J.A., 2003. A generalized single-channel method for retrieving land surface temperature from remote sensing data. J. Geophys. Res.: Atmos. 108. Jönsson, P., Eklundh, L., 2004. Timesat — a program for analyzing time-series of satellite sensor data. Comput. Geosci. 30, 833–845. Kljun, N., Calanca, P., Rotach, M., Schmid, H., 2015. A simple two-dimensional parameterisation for flux footprint prediction (ffp). Geosci. Model Dev. 8, 3695–3713. Kustas, W.P., Diak, G.R., Norman, J.M., 2001. Time difference methods for monitoring regional scale heat fluxes with remote sensing. Land Surface Hydrology, Meteorology, and Climate: Observations and Modeling. pp. 15–29. Kustas, W.P., Norman, J.M., 1999. Evaluation of soil and vegetation heat flux predictions using a simple two-source model with radiometric temperatures for partial canopy cover. Agric. For. Meteorol. 94, 13–29. Kustas, W.P., Norman, J.M., 2000a. Evaluating the effects of subpixel heterogeneity on pixel average fluxes. Remote Sens. Environ. 74, 327–342. Kustas, W.P., Norman, J.M., 2000b. A two-source energy balance approach using directional radiometric temperature observations for sparse canopy covered surfaces. Agron. J. 92, 847–854. LSA-SAF, 2011. LSA-SAF Validation Report Down-welling Surface Shortwave Flux (DSSF). Products: LSA-07 (MDSSF), LSA-08 (EDSSF), LSA-09 (DIDSSF). Moore, C., 1986. Frequency response corrections for eddy correlation systems. BoundaryLayer Meteorology 37, 17–35. Norman, J.M., Kustas, W.P., Humes, K.S., 1995. Source approach for estimating soil and vegetation energy fluxes in observations of directional radiometric surface temperature. Agric. For. Meteorol. 77, 263–293. Pasolli, L., Asam, S., Castelli, M., Bruzzone, L., Wohlfahrt, G., Zebisch, M., Notarnicola, C., 2015a. Retrieval of leaf area index in mountain grasslands in the alps from modis satellite imagery. Remote Sens. Environ. 165, 159–174. Pasolli, L., Notarnicola, C., Bertoldi, G., Bruzzone, L., Remelgado, R., Greifeneder, F., Niedrist, G., Della Chiesa, S., Tappeiner, U., Zebisch, M., 2015b. Estimation of soil moisture in mountain areas using SVR technique applied to multiscale active radar images at c-band. IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens. 8, 262–283. Priestley, C., Taylor, R., 1972. On the assessment of surface heat flux and evaporation using large-scale parameters. Mon. Wea. Rev 100, 81–92. Santanello, J.A., Friedl, M.A., 2003. Diurnal covariation in soil heat flux and net radiation. J. Appl. Meteorol. 42, 851–862. Schotanus, P., Nieuwstadt, F., De Bruin, H., 1983. Temperature measurement with a sonic anemometer and its application to heat and moisture fluxes. Bound.-Lay. Meteorol. 26, 81–93. Sobrino, J.A., Jiménez-Muñoz, J.C., Paolini, L., 2004. Land surface temperature retrieval from landsat TM 5. Remote Sens. Environ. 90, 434–440. Stöckli, R., 2013. The HelioMont surface solar radiation processing. Bundesamt für Meteorologie und Klimatologie, MeteoSchweiz. Stoy, P.C., Mauder, M., Foken, T., Marcolla, B., Boegh, E., Ibrom, A., Arain, M.A., Arneth, A., Aurela, M., Bernhofer, C., et al., 2013. A data-driven analysis of energy balance closure across fluxnet research sites: the role of landscape scale heterogeneity. Agric. For. Meteorol. 171, 137–152. Twine, T.E., Kustas, W., Norman, J., Cook, D., Houser, P., Meyers, T., Prueger, J., Starks, P., Wesely, M., 2000. Correcting eddy-covariance flux underestimates over a grassland. Agric. For. Meteorol. 103, 279–300. Wan, Z., Zhang, Y., Zhang, Q., Li, Z.-L., 2004. Quality assessment and validation of the modis global land surface temperature. Int. J. Remote Sens. 25, 261–274. Webb, E.K., Pearman, G.I., Leuning, R., 1980. Correction of flux measurements for density effects due to heat and water vapour transfer. Q. J. R. Meteorol. Soc. 106, 85–100. Wilczak, J.M., Oncley, S.P., Stage, S.A., 2001. Sonic anemometer tilt correction algorithms. Bound.-Lay. Meteorol. 99, 127–150. Wohlfahrt, G., Albin, H., Georg, N., Katharina, S., Enrico, T., Peng, Z., 2016. On the energy balance closure and net radiation in complex terrain. Agric. For. Meteorol. 226, 37–49. Wohlfahrt, G., Anfang, C., Bahn, M., Haslwanter, A., Newesely, C., Schmitt, M., Drösler, M., Pfadenhauer, J., Cernusca, A., 2005. Quantifying nighttime ecosystem respiration of a meadow using eddy covariance, chambers and modelling. Agric. For. Meteorol. 128, 141–162. Wohlfahrt, G., Haslwanter, A., Hörtnagl, L., Jasoni, R.L., Fenstermaker, L.F., Arnone, J.A., Hammerle, A., 2009. On the consequences of the energy imbalance for calculating surface conductance to water vapour. Agric. For. Meteorol. 149, 1556–1559. Wohlfahrt, G., Irschick, C., Thalinger, B., Hörtnagl, L., Obojes, N., Hammerle, A., 2010. Insights from independent evapotranspiration estimates for closing the energy balance: a grassland case study. Vadose Zone J. 9, 1025–1033. Wohlfahrt, G., Widmoser, P., 2013. Can an energy balance model provide additional constraints on how to close the energy imbalance? Agric. For. Meteorol. 169, 85–91.

Martin, P., Berbigier, P., Bernhofer, C., et al., 2000. Estimates of the annual net carbon and water exchange of forests: the EUROFLUX methodology. vol. 30. Baldocchi, D.D., Hincks, B.B., Meyers, T.P., 1988. Measuring biosphere-atmosphere exchanges of biologically related gases with micrometeorological methods. Ecology 1331–1340. Bastiaanssen, W., Menenti, M., Feddes, R., Holtslag, A., 1998. A remote sensing surface energy balance algorithm for land (sebal). 1. formulation. J. Hydrol. 212, 198–212. Bertoldi, G., Chiesa, S.D., Notarnicola, C., Pasolli, L., Niedrist, G., Tappeiner, U., 2014. Estimation of soil moisture patterns in mountain grasslands by means of {SAR} {RADARSAT2} images and hydrological modeling. J. Hydrol. 516, 245–257. http:// www.sciencedirect.com/science/article/pii/S0022169414001255http://dx.doi.org/ 10.1016/j.jhydrol.2014.02.018. Determination of soil moisture: measurements and theoretical approaches. Brutsaert, W., Sugita, M., 1992. Application of self-preservation in the diurnal evolution of the surface energy budget to determine daily evaporation. J. Geophys. Res.: Atmos. (1984-2012) 97, 18377–18382. Burba, G.G., McDERMITT, D.K., Grelle, A., Anderson, D.J., Xu, L., 2008. Addressing the influence of instrument surface heat exchange on the measurements of CO2 flux from open-path gas analyzers. Glob. Chang. Biol. 14 (8), 1854–1876. Cammalleri, C., Anderson, M., Gao, F., Hain, C., Kustas, W., 2014a. Mapping daily evapotranspiration at field scales over rainfed and irrigated agricultural areas using remote sensing data fusion. Agric. For. Meteorol. 186 (0), 1–11. http://www. sciencedirect.com/science/article/pii/S0168192313002815http://dx.doi.org/10. 1016/j.agrformet.2013.11.001. Cammalleri, C., Anderson, M., Kustas, W., 2014b. Upscaling of evapotranspiration fluxes from instantaneous to daytime scales for thermal remote sensing applications. Hydrol. Earth Syst. Sci. 18, 1885–1894. Castelli, M., Stöckli, R., Zardi, D., Tetzlaff, A., Wagner, J., Belluardo, G., Zebisch, M., Petitta, M., 2014. The heliomont method for assessing solar irradiance over complex terrain: validation and improvements. Remote Sens. Environ. 152, 603–613. Fang, H., Wei, S., Liang, S., 2012. Validation of modis and cyclopes lai products using global field measurement data. Remote Sens. Environ. 119, 43–54. Fisher, J.B., Tu, K.P., Baldocchi, D.D., 2008. Global estimates of the land-atmosphere water flux based on monthly AVHRR and ISLSCP-II data, validated at 16 fluxnet sites. Remote Sens. Environ. 112, 901–919. Foken, T., 2008. The energy balance closure problem: an overview. Ecol. Appl. 18 (6), 1351–1367. Gage, S.H., Doll, J.E., Safir, G.R., 2015. A Crop Stress Index to Predict Climatic Effects on Row-crop Agriculture in the US North Central Region. The Ecology of Agricultural Landscapes: Long-Term Research on the Path to Sustainability. pp. 77–103. Gan, G., Gao, Y., 2015. Estimating time series of land surface energy fluxes using optimized two source energy balance schemes: model formulation, calibration, and validation. Agric. For. Meteorol. 208, 62–75. Gao, F., Anderson, M.C., Kustas, W.P., Wang, Y., 2012. Simple method for retrieving leaf area index from landsat using modis leaf area index products as reference. J. Appl. Remote. Sens. 6 063554–1. Gobiet, A., Kotlarski, S., Beniston, M., Heinrich, G., Rajczak, J., Stoffel, M., 2014. 21st century climate change in the European Alps — a review. Sci. Total Environ. 493, 1138–1151. http://www.sciencedirect.com/science/article/pii/ S0048969713008188http://dx.doi.org/10.1016/j.scitotenv.2013.07.050. Guillevic, P.C., Hulley, G.C., Hook, S.J., Olioso, A., Sanchez, J.M., Drewry, D., Running, S.W., Fisher, J.B., 2014. Evapotranspiration from Airborne Simulators as a Proxy Datasets for NASA's ECOSTRESS Mission — A New Thermal Infrared Instrument on the International Space Station. AGU Fall Meeting Abstracts. 1. pp. 0802. http:// adsabs.harvard.edu/abs/2014AGUFM.H41B0802G. Guzinski, R., Anderson, M., Kustas, W., Nieto, H., Sandholt, I., 2013. Using a thermalbased two source energy balance model with time-differencing to estimate surface energy fluxes with day-night modis observations. Hydrol. Earth Syst. Sci. 17, 2809–2825. Hain, C.R., Crow, W.T., Anderson, M.C., Mecikalski, J.R., 2012. An ensemble Kalman filter dual assimilation of thermal infrared and microwave satellite observations of soil moisture into the Noah land surface model. Water Resour. Res. 48. Hammerle, A., Haslwanter, A., Schmitt, M., Bahn, M., Tappeiner, U., Cernusca, A., Wohlfahrt, G., 2007. Eddy covariance measurements of carbon dioxide, latent and sensible energy fluxes above a meadow on a mountain slope. Bound.-Layer Meteorol. 122, 397–416. Haslinger, K., Schöner, W., Anders, I., 2015. Future drought probabilities in the greater alpine region based on cosmo-clm experiments — spatial patterns and driving forces. Meteorologische Zeitschrift 25 (2), 137–148. Haslwanter, A., Hammerle, A., Wohlfahrt, G., 2009. Open-path vs. closed-path eddy covariance measurements of the net ecosystem carbon dioxide and water vapour exchange: a long-term perspective. Agric. For. Meteorol. 149, 291–302. Hiller, R., Zeeman, M.J., Eugster, W., 2008. Eddy-covariance flux measurements in the complex terrain of an alpine valley in Switzerland. Bound.-Lay. Meteorol. 127, 449–467. Hochberg, E.J., Roberts, D.A., Dennison, P.E., Hulley, G.C., 2015. Special issue on the hyperspectral infrared imager (hyspiri): emerging science in terrestrial and aquatic ecology, radiation balance and hazards. Remote Sens. Environ. 167, 1–5. http:// www.sciencedirect.com/science/article/pii/S0034425715300420http://dx.doi.org/ 10.1016/j.rse.2015.06.011. Special Issue on the Hyperspectral Infrared Imager

342