Using uncertainty of Penman and Penman–Monteith methods in combined satellite and ground-based evapotranspiration estimates

Remote Sensing of Environment 169 (2015) 102–112

Contents lists available at ScienceDirect

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

Using uncertainty of Penman and Penman–Monteith methods in combined satellite and ground-based evapotranspiration estimates R.S. Westerhoff Deltares, Delft, The Netherlands University of Waikato, Hamilton, New Zealand GNS Science, Taupo, New Zealand

a r t i c l e

i n f o

Article history: Received 15 March 2015 Received in revised form 12 June 2015 Accepted 16 July 2015 Available online 21 August 2015 Keywords: Evapotranspiration MOD16 Ground-based estimates Uncertainty Sensitivity analysis Monte-Carlo fitting Reference ET PET

a b s t r a c t Satellite data are often used for their ability to fill in temporal and spatial patterns in data-sparse regions. It is also known that global satellite products generally contain more noise than ground-based estimates. Data validation of satellite data often treats ground-based estimates as the ‘gold standard’: without error or uncertainty. In the estimation of evapotranspiration (ET) however, ground-based estimates have considerable uncertainty, caused by the input components of the ET equations. This research presents an analysis of uncertainty of reference ET (ET0) caused by these input components. A dataset of correlated random variables is generated for a country with a diverse climate and diverse density of ground observations: New Zealand. The uncertainty analysis shows that: ET0 is most sensitive to temperature, followed by solar radiation, relative humidity, and cloudiness ratio; and that uncertainty varies between 10% and 40% of ET0, and depends on the ET0 value. Using this uncertainty analysis, a set of correlated random variables, and a Monte-Carlo fitting approach, MOD16 satellite PET data becomes a ‘soft interpolator’ between ground-based ET0 estimates. The resulting 1 km × 1 km monthly nation-wide dataset has the advantage of: taking into account land cover and vegetation characteristics through the use of satellite data; still abiding to local climate diversity and locally used standards through the use of ground-based estimates; and containing an uncertainty estimate. Further comparison suggests that original MOD16 satellite PET could estimate real PET better than using ground-based estimates of ET0. Further research recommends combination with other existing gridded ET estimates, and further validation of real PET estimates. © 2015 The Author. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction Evapotranspiration (ET) can be defined as liquid water from soil transferring to vapour in the atmosphere (Irmak, 2012). It can either be expressed as potential ET (PET: the amount of ET if ample water is available) or actual ET (AET). ET is an important yet complex part of the water cycle, as it cannot be measured directly, but is inferred through energy budget methods (e.g., Anderson, Norman, Mecikalski, Otkin, & Kustas, 2007; Bastiaanssen et al., 1998) or hydrological water balance estimations, with measurements of vegetation status and climate variables (i.e., stomatal conductance, temperature, solar radiation, air pressure, air humidity, wind; Mu, Zhao, & Running, 2011). PET is an important parameter, as hydrologists often estimate PET first before further deriving the more complex AET. More often, PET is first simplified to a reference ET (ET0), which provides a PET estimate for a standard condition (e.g., green grass of 12 cm, Allen, 1998). ET0 estimates are generally ‘easier’ to estimate using observed meteorological variables at climate stations, compared to the more expensive and work-intensive ET estimates through lysimeters or flux towers. ET can be a large yet diverse part of the water budget. For example, in New Zealand, AET is on average approximately 53% of mean

annual rainfall, but varies considerably throughout the country: mean annual values range from 100 to 1100 mm/year (Woods, Hendrikx, Henderson, & Tait, 2006). Knowing the uncertainty of AET, PET and/or ET0 is important for hydrologists, as it is such a large part of the water budget. For example, Howard and Lloyd (1979) showed that 10% overestimation of AET leads to 5% underestimation of yearly groundwater recharge and up to 20% in summer groundwater recharge. This research does not provide a description of the many variations of ET estimation methods, as that is already given in the comprehensive overview of McMahon, Peel, Lowe, Srikanthan, and McVicar (2013). In this study, ET0 estimates are inferred through Penman (1963) and Monteith (1965); Penman (1948) methods. 1.1. Global satellite products could benefit ET estimation on a catchment scale In recent years, global satellite data products are increasingly developed. Space research on water resources, drought and floods has been increasingly linked into collaborative global efforts to develop data on a global scale. Examples are large-scale efforts and projects such as GEWEX (Global Energy and Water Exchange; GEWEX, 2014), and GEO

http://dx.doi.org/10.1016/j.rse.2015.07.021 0034-4257/© 2015 The Author. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

(Group on Earth Observation; GEO, 2014b), who oversee and coordinate global climate and space research for water (GEO, 2014a). Recent research projects of the European 7th Framework Programme (i.e., eartH2Observe, 2015; GLOWASIS, 2015) aim at providing open data on global water resources datasets for scientific communities. These communities could then use these data for research on their continent, country or catchment. As most policy on water is on a national and regional scale, global satellite data could thus be adding value in countries where ground observations are sparse. However, most global satellite data and derived products are not easily projected to a country or catchment scale. For example, for small catchments, some global ET data (e.g., from Miralles, De Jeu, Gash, Holmes, & Dolman, 2011; Miralles et al., 2011) is often too coarse. Also, most satellite-derived ET products embed a ‘pseudo-model’, which can hold over-simplified model assumptions and unknown uncertainty (e.g., Miralles, De Jeu et al., 2011; Miralles, Holmes et al., 2011; Mu, Heinsch, Zhao, & Running, 2007; Mu et al., 2011). The oneon-one use of these data on a catchment scale is thus prone to errors. A future challenge lies ahead: how to use the large amount of available global satellite ET data on a catchment scale without the results becoming unreliable. 1.2. Treating ground-based estimates as the ‘gold standard’ is a wrong assumption Simple interpolation methods would lead to a gridded product if ample ground-based data are available, the ground-based value logically being treated as the correct value. However, this is not always the case, as ground-based networks are usually sparse and can furthermore also have uncertainty. In the case of sparse data, additional (e.g., satellite) data can be used. This can be done with kriging techniques, e.g., described in Chirlin and Wood (1982), who proposed universal kriging for merging or assimilating in-situ station data into a gridded product. Interpolation methods that use other variables as an additional constraint are also used in data-sparse environments. For example, Tait and Woods (2007) use elevation as an additional parameter and describe 5 km × 5 km daily ET0 for New Zealand by interpolating long term daily data of 70 climate stations using trivariate thinsplicing interpolation. These techniques correlate well with groundbased estimates, but they do not take into account the uncertainty of those ground-based ET estimates. In most studies estimating ET from satellite data (Bastiaanssen et al., 1998; Miralles, Holmes, et al., 2011; Mu et al., 2007, 2011; Su, 2002) the satellite estimations are validated with ground-based estimates from (micro-) meteorological or eddycovariance flux towers. Validations seldom lead to a one-on-one correlation, and it is considered in these studies that the ground-based estimates are correct and that the satellite estimates contain noise and uncertainty. Earth-observing satellites measure from several hundreds of km above where the atmosphere creates incoherent noise. It is also known that satellite ET data can have high errors in areas with high elevation differences (Van der Tol & Parodi, 2012). However, groundbased estimates are also uncertain: it is for example known that flux tower ET estimations can lead to underestimation of the energy balance by 20% (Wilson et al., 2002). Also, most PET estimations from climate stations use simplifying assumptions. For example, FAO56 ET0 by Allen (1998) assumes grass with a constant height of 0.12 m, which for a well-watered crop different from that grass should then be corrected by a crop coefficient to get the real PET for that crop: crop evapotranspiration. Errors due to different vegetation types or wrong vegetation growth status are often not considered, and the effect on real PET is relatively unknown, but potentially large. Satellite data are able to measure vegetation characteristics and could thus help in a better PET estimation. Ideally, one wants a combination of both satellite ET and ground-based estimates, which takes uncertainty of both datasets into account. Potentially, weak points from both methods could then be scored out.

103

1.3. It is hard to find a uniform uncertainty or sensitivity analysis that takes into account time-scales, seasonality and data inter-dependency Uncertainty comes in different forms and is therefore difficult to assess on a uniform spatial and temporal scale. It could stem from measurements, or models that ignore inter-dependency of input components. Literature values of uncertainties cannot always be used in further studies, as they: only represent the local or regional climate and observation strategies; they comprise different time windows (e.g., daily, monthly or annual); have different model origin; or have different statistical origin. Error and sensitivity analyses of ET to its input components have been done in the past, but with different outcomes. For example, Saxton (1975) showed that computed ET was most sensitive to changes or errors in net radiation: in summer, 50–90% was governed by net radiation, and 20–30% by aerodynamic properties (vapour pressure deficit, wind). Gong, Xu, Chen, Halldin, and Chen (2006) showed that ET was most sensitive to relative humidity, followed by solar (shortwave) radiation, air temperature and wind speed and that sensitivities changed over the season. Coleman and DeCoursey (1976) estimated standard deviation of several ET methods and found average errors of approximately 0.2 mm/day in winter and 3 mm/day in summer. However, these were for different climatic environments and different ET estimation methods. Estévez, Gavilàn, and Berengena (2009) show that in Southern Spain ET is most sensitive to temperature in summer and it causes overestimation of ET; relative humidity in winter causes underestimation of ET. All above-mentioned studies show that sensitivity of ET to its input components depends on climate region and season. Uncertainty of some ET input components has been described by Tait and Zheng (2007). They established a range of model prediction errors of minimum and maximum median annual temperature of 0.4 and 0.8 °C using the ANUSPLIN method (Australian National University, 2014). For relative humidity, Lin and Hubbard (2004) found uncertainties, both modelled and measured, to be in the range of 0.5 to 5%, depending on the temperature. For wind speed Leathwick and Stephens (1998) used two different methods and found uncertainties in between 0.6 and 1.6 m/s for daily wind speed measurements; Tait and Zheng (2007) found uncertainties in between 0.6 and 1.1 m/s for median annual wind speed. Solar radiation measured from earth's surface is commonly used for measuring cloudiness or surface albedo, but knows uncertainty of site measurements and down-scaling methods Li, Cribb, and Chang (2005). Surface albedo values are derived from solar radiation measurements, and therefore complex, because these again are correlated with cloudiness. The inter-relation of these parameters implies that they are not independent. There is a need for better uncertainty assessment of ET and its many components, taking into account this inter-dependence. This study first describes a method for estimating the uncertainty of ground-based ET0 estimates and its different input components. This uncertainty is then used to project a global satellite PET product to ET0, abiding the locally used estimates and methods. The methods and their results are explained for monthly historical ET0 in a case study in New Zealand. Direct application of original satellite PET using the full P–M method is furthermore discussed. 2. Input data 2.1. Satellite PET The MOD16 algorithm uses satellite data from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) sensor onboard the Terra and Aqua satellites. The satellite-derived parameters are land cover, albedo, leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FPAR), and enhanced vegetation index (EVI). Temperature, radiation, air humidity and pressure data are derived from daily global meteorological reanalysis data set from NASA's Global Modeling and Assimilation Office (GMAO). The algorithm, described in

104

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

Mu et al. (2007) and Mu et al. (2011), uses the P–M approach. Data are available online in HDF (Hierarchical Data Format) files (NTSG, 2013). The spatial resolution of the cells is 0.926 km by 0.926 km, but was re-gridded to 1 km × 1 km for this study. Data are available in 8-daily, monthly, and yearly intervals. The monthly data were used in this study. Four variables are available in the HDF file: PET; AET; potential latent heat flux; and latent heat flux. Mean annual MOD16 PET, compiled from monthly data, is shown in Fig. 2. MOD16 data currently covers the period January 2000 to December 2014. For this study monthly MOD16 PET was used as input in the method, as it could be best evaluated against groundbased estimates. These data will be called original MOD16 PET onwards. 2.2. Ground-based ET0 Reference ET ground-based estimates from 112 climate stations (Fig. 1) were downloaded from the New Zealand's Climate Database web portal CliFlo (NIWA, 2014a,b). Data were compiled to monthly values from January 2000 to December 2014. Some stations do not cover 2000–2014 entirely, and there are also gaps for some years and months in the data. Data from the climate stations have been processed using a Penman reference ET method based on Penman (1963). Conversion to the more widely used FAO56 reference ET of Allen (1998) is currently planned. Both methods are described in more detail in Section 3. 2.3. Meteorology data Daily climate data from 2000 to 2014 for New Zealand are available in a regular ~5 km grid from the Virtual Climate Station (VCS) network (NIWA, 2014a,b; Tait, 2014). The VCS data used in this research were rainfall according to Tait, Henderson, Turner, and Zheng (2006); relative humidity; solar radiation; maximum and minimum temperature; and wind speed. 3. Theory and methods 3.1. Estimating uncertainty of ground-based ET0 estimates 3.1.1. Theory of Penman–Monteith ET estimation The Penman–Monteith (P–M) equation (Monteith, 1965; Penman, 1948) calculates ET based on the earth's surface energy balance and

Fig. 2. PET from MOD16 data. Mean annual value for 2000–2014.

atmospheric water demand. Commonly used simplified equation to calculate reference ET (Allen, 1998; Penman, 1963) calculate PET by using simplified assumptions of land cover and growth. Reference ET according to Allen (1998) is most used and is also known as ‘FAO56 reference ET’ or ‘FAO56 reference crop ET’. It is called ET0,FAO56 onwards. The Penman method, based on the same input components, is mostly used in New Zealand as an alternative to the ET0,FAO56. It calculates evaporation from a grass surface (Burman & Pochop, 1994; NIWA, 2014a,b). Although not accounting for transpiration, Penman ET0 is used as such, and will be called ET0,Penman onwards. Both ET0,FAO56 and ET0,Penman, can be simplified to a function of limited input components: 60:3 u2 ðes −ea Þ T þ 273 ¼ 0:408  Δ þ 0:067ð1 þ 0:34u2 Þ ΔRn þ

ET 0; FAO56

ð1Þ

and ET 0;Penman ¼

ΔRn þ 0:4244ð1 þ 0:536u2 Þðes −ea Þ ðΔ þ 0:067Þð2:501−0:00236T Þ

ð2Þ

where ET is in [mm day−1]; u2 is wind speed at 2 m above the surface [m s−1]; T is temperature [°C]; es is saturation vapour pressure at the surface [kPa]; ea is actual vapour pressure of the air above the evaporating surface [kPa]; – Δ is the gradient of the saturation pressure curve [kPa °C−1]; – RH is the relative humidity [%]; – Rn is the net radiation of the earth's surface [MJ m−2 day−1].

– – – – –

Fig. 1. Location of 112 meteorological stations, which were used to estimate ET0 in this study.

Eq. (1) simplifies transpiration for a green grass with a height of 0.12 m, while Eq. (2) does not take into account transpiration of the plant through stomata. Both equations assume no soil water deficit. Several approaches exist for calculating the input components. The approach used here is explained in more detail in the Annex. For this study, T, RH, S0, u2, n/N and α are called core input components, as Rn, Δ, es and ea are generally derived from them.

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112 Table 1 Values of components in the FAO56 and Penman equations through literature and multivariate distributions of VCS data. p1 = 1st percentile, p99 = 99th percentile. p1

p99

Units

23.3 96.2 50.7 6.26 0.40 0.82

°C % MJ m−2 day−1 m s−1 – –

Estimated through core input es 0.75 2.86 Δ 0.05 0.17 Sn 1.03 25.4 Ln 1.10 2.95 Rn 0 22.8

kPa kPa °C−1 MJ m−2 day−1 MJ m−2 day−1 MJ m−2 day−1

Core input T 2.92 RH 69.8 S0 3.70 u2 0.87 α 0.08 n/N 0.07

Reference

Hendriks (2010) Tait and Zheng (2007) Ahrens (2007); Oke (1992)

3.1.2. Method of estimation of ET0 uncertainty For both ET0,Penman and ET0,FAO56, a set of one million correlated random realisations of the core input components were generated. Correlated random realisations of d components fall within a multivariate normal distribution f (Mathworks, 2014a; Rose & Smith, 1996): −1 1 1 f ðx; μ; VÞ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e2ðx−μ ÞV ðx−μ Þ ; d jVjð2π Þ

ð3Þ

with f a function of 1-by-d vectors x (locations within distribution), μ (means of input components), and d-by-d variance-covariance matrix V. The input components were T, RH, S0, u2, α and n/N. Before applying Eq. (3), these pre-processing steps were taken: – Daily VCS Tmin, Tmax, RH, u2 and solar radiation Rad were compiled to monthly values; – T was estimated as T = (Tmin + 2Tmax)/3 – Monthly S0 throughout the year was estimated from Hendriks (2010) (his Figure B2.12.2); – Although albedo values for the FAO56 and Penman methods are fixed (0.23 and 0.25, respectively), they were randomly distributed in between 0.08 (water) and 0.4 (dry sand), to look at the sensitivity of the methods to this input component; – n/N was calculated using Eqs. (A.6) and (A.7): n 2Rad ¼ −0:5: N S0

ð4Þ

– μ and V were calculated.

Table 2 Estimated uncertainties in input components of FAO56 and Penman equations. p1 = 1st percentile, p99 = 99th percentile. p99

Units

Core input components σT 0.01 σRH 0.55 σ S0 0 σ u2 0.51 σα 0 σn/N 0.05

p1

0.99 4.96 0 1.49 0.01 0.15

°C % MJ m−2 day−1 m s−1 – –

Propagated with core input σ es 0 σΔ 0 σ Sn 0.15 σ Ln 0.30 σ Rn 0

0.14 0.01 3.18 0.46 3.20

kPa kPa °C−1 MJ m−2 day−1 MJ m−2 day−1 MJ m−2 day−1

105

Table 3 Sample sizes, slopes A, intercepts B and RMSE for all monthly linear fits.

Average Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec

Sample size

A

σA

B

σB

1087 1065 1078 1064 1082 1087 1085 1093 1105 1110 1092 1095 1083

1.02 0.88 1.12 1.31 1.28 1.06 0.93 0.87 1.05 1.25 1.02 0.80 0.70

0.08 0.09 0.09 0.08 0.07 0.05 0.05 0.05 0.07 0.12 0.11 0.11 0.11

27.8 66.6 28.2 4.4 6.4 11.7 13.0 14.0 9.8 2.47 27.0 61.7 88.2

6.7 11.9 9.1 7.3 3.3 1.3 0.8 1.0 2.3 6.8 9.2 13.2 14.6

Values and uncertainties for further input components in es, ea, Δ, Rn were then calculated with Eqs. (A.1) to (A.5). This was done assuming basic error propagation (Tellinghuisen, 2001): σ f 2 ¼ gT Vg

ð5Þ

where σ f 2 is the variance of a function f, which has n = 1:N input components; g is a vector of input component ∂f/∂ni; and gT is the transpose of g. All resulting input components were then used to generate one million realisations of ET0. Then, a million random realisations of uncertainties in the core input components (from the variances on the diagonal of V, according to Table 1) were used to calculate a ‘noisy’ ET0. The difference between the normal and noisy ET0 is called the error in ET0 due to propagated error of the input component, or ϵET0. The standard deviation of the mean values is called the uncertainty, or σET0. All realisations of σET0 were then normalized to ET0. The ET 0–σET0 relation was then averaged per ET 0 value (bin size ~ 0.33) and fitted to a function with a smoothing spline fit (Mathworks, 2014b). This function, describing the uncertainty of in-situ ET, was used in the further calculations. For further processing, it was assumed that the maximum uncertainty of satellite PET is larger than that of ground-based estimates, as was already explained in the introduction. It is assumed that satellite uncertainty is 150% of groundbased uncertainty, as there is no other value known for such as en error estimate (more in discussion). 3.2. Projection of satellite PET to ET0 ground-based estimates A linear relation between original MOD16 PET and ground-observed ET0 data was calculated for monthly data using a Monte-Carlo approach. This was done by least-squares fitting all available original MOD16 PET with ground-based ET0 from all available meteorological stations using 1000 realisations of error in both datasets. Each realisation used a random deviation within the uncertainties of both datasets. When the linear fit is Table 4 Specifications of MOD16 ET data files. Original MOD16 data has been regridded. Format Spatial resolution Coordinate system Start End Time interval Unit Data layers

NetCDF UCAR (2014), CF-1.6 30 arc sec WGS84 (EPSG:4326) January 2010 December 2014 Monthly mm/month Monthly ET MOD160 Uncertainty in MOD160 Original MOD16 PET Original MOD16 AET

106

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

Fig. 3. Means and standard deviation σ for ET0,FAO for correlated random values of the ET core input components T, RH, u2, n/N, albedo α and S0. εnoise is the error when noise in the input components is introduced.

defined as y = Ax + B, the realisation with the lowest χ2 coefficient was chosen as the best fit. χ2 is defined by Taylor (1997) as Eq. (6): χ2 ¼

N X ðy −B−Axi Þ2 i

i¼1

σ 2y

;

ð7Þ

in which σA and σB are the standard deviations in slope and intercept, respectively. These were derived by calculating the standard deviation of the χ2 distribution. To account for seasonal variation, original MOD16 PET for each month was fitted to ground-observed ET0. This led to 12 solutions (January to December) for slope and intercept. These were then used to project original MOD16 PET to ground-observed ET0 using slope and intercept of the best linear fit for each month: MOD160 ¼

MOD16orig −B ; A

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σ MOD16orig 2 þ σ A 2 þ σ B 2 :

ð9Þ

ð6Þ

with: x representing the ground-based estimates; y the original MOD16 PET data; σy the total uncertainty for y; and N the number of groundbased estimates. Using the A and B of the best fit, the original MOD16 PET was then defined as Eq. (7): MOD16orig ¼ ðA  σ A ÞET 0 þ B  σ B ;

with estimated uncertainty:

ð8Þ

The resulting data were gridded to an openly available gridded NetCDF file (UCAR, 2014, further specifications are shown in Table 4). The original MOD16 PET and AET data were regridded in the same format. 4. Results 4.1. Estimation of ET0 uncertainty Values and uncertainties for all ET input components, resulting from the multivariate random distributions, are shown as percentiles 1 and 99 (p1 and p99) in Tables 1 and 2. Fig. 3 shows that values of ET0,FAO range from approximately 0 to 10 mm/day (for Penman 0 to 15 mm/day, Fig. 4), and that ET0 is most sensitive to its input components temperature, followed by cloudiness ratio and relative humidity (the slope of the red line). The grey band denotes the standard deviation of ET0 outcomes of all random generations, while the black-dotted line gives an estimate of the effect of noise per input component. Some plots show a ‘saw-tooth’ like pattern (e.g. at high RH in Fig. 4, top

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

107

Fig. 4. Means and standard deviation σ for ET0,Penman for correlated random values of the ET core input components T, RH, u2, n/N, albedo α and S0. εnoise is the error when noise in the input components is introduced.

right), which is caused by the less available values in the data bin. The above results can be used in further sensitivity analyses. However, the result used for further processing is the mean error resulting from all

noise measurements per value of ET0 (denoted as σET0 onwards, and shown for FAO56 and Penman ET0 in Figs. 5 and 6, respectively). The ratio of σET0 to ET0 approaches 10% for high ET0, and 40% for low ET.

Fig. 5. Means and standard deviation σ of the normalized uncertainty σET0 (FAO56) for correlated random values of T, RH, u2, n/N, albedo α and S0.

Fig. 6. Means and standard deviation σ of the normalized uncertainty σET0 (Penman) for correlated random values of T, RH, u2, n/N, albedo α and S0.

108

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

Penman results show that ET0 is more sensitive and uncertain for high ET0 values, reaching ratios of almost 30%. 4.2. Projection of satellite PET to ET0 ground-based estimates Original MOD16 PET are plotted against ground-based estimates for four months (January, April, July, Oct, the ‘middle of the season months’) in 2010 in Fig. 7. Original MOD16 PET is generally higher than groundobserved data. This is as expected: the P–M approach takes into account vegetation transpiration, where the Penman approach does not. The monthly Monte-Carlo least squares fit for each month (Fig. 8 and Table 3) show a seasonal pattern, where warmer months show different fits than colder months. The projected MOD160 data fit well with most ground-based estimates (Fig. 9). Comparison of MOD160 data for random years at random ground-based estimates, shows that MOD160 data, although comparing well, can show some distinct seasonal differences (Fig. 10). This is also as expected: MOD16 data embeds land cover data with different vegetation. The resulting MOD160 ET0 dataset, as well as its uncertainty, is a monthly 1 km × 1 km product. A compiled image of the data shows mean annual ET0 over the period 2000–2014 (Fig. 11), including its uncertainty (Fig. 12).

5. Discussion The main advantage of the presented monthly MOD160 ET0 dataset is that it has a high resolution (1 km × 1 km), has uncertainty information, and takes into account vegetation characteristics (as measured by MODIS), while still abiding to the locally used standards in groundbased estimates. Other gridded ET estimates (i.e., 5 km × 5 km daily ET estimates from Tait & Woods, 2007) are interpolated based on elevation (see Section 1). Elevation is indeed related to ET input component temperature and solar radiation, but less obvious to e.g. albedo, wind speed or land cover. A potential drawback of MOD16 satellite data is that input data on global meteorology are from global sources and could maybe not capture climate diversity of countries like New Zealand, making the data less reliable. This was already concluded by Zhao, Running, and Nemani (2006), who evaluated the sensitivity of MODIS global terrestrial gross and net primary (carbon) production to the uncertainties of meteorological inputs on global vegetated areas and found that biases in meteorological reanalyses can introduce substantial error to these estimations. However, the comparison with original MOD16 data and Penman ground-based estimates in colder months (Fig. 7) shows that original MOD16 data is just slightly higher than Penman ET. This is as expected: a slightly higher ET in relatively colder months shows the minor

Fig. 7. Original MOD16 PET for January, April, July and October of 2010 compared to ground-based ET0 estimates. The stations are plotted from the southernmost (left) to northernmost (right) locations.

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

109

Fig. 8. Monthly best fits for original MOD16 PET to ground-based ET0 estimates. Black dashed line: best fit per month. Grey lines: linear fits within the standard deviation of slope and intercept. Error bars are not shown.

difference that the transpiration causes. The qualitative trend shows this difference increasing in warmer months, which seems to be correct. Although having a potentially larger error, the original MOD16 PET data could be used as a first estimate of PET, as it is the only estimate taking into account vegetation characteristics. This could potentially enhance the capability of this dataset to be used for irrigation studies (i.e., the data could be used to analyse crop growth in irrigated areas). However, there are no ground-based estimates through which the original MOD16 PET can be validated, except for the mentioned correct qualitative trend. Satellite MOD16 data as presented in this research could be used to improve existing daily (VCS) estimates: merge of MOD160 ET0 datasets with those from Tait and Woods (2007), followed by a second merge to PET estimates using the original MOD16 data, could help score out the potential drawbacks of each method, and could aid in better estimation of crop coefficient factors using the best-fit values from Table 3. A conversion of ET ground-based estimates from Penman to FAO56 would help improve this merge. The FAO56, Penman, and MOD160 ET0 estimates do not take into account soil moisture deficit (as they are reference ET). Soil moisture deficit is represented in the P–M equation in the parameter surface resistance rs, which has complex inter-relations with soil moisture, and stomatal conductance of different plant and tree types (Verhoef & Egea, 2014). As ground-based estimates calculate a reference ET, this complexity was not tackled in this research. Suggestions for a simplified conversion of ET0 to AET are: (1) to use the original MOD16 AET; (2) to use the original MOD16 PET in a hydrological soil water balance approach (e.g., Thornthwaite, 1948; Thornthwaite & Mather, 1957) to estimate AET based on soil water deficit from the soil storage component. As satellite data in this study were processed to ET0, the error analysis is therefore also simplified. Recommended further consideration of ET uncertainty analyses could therefore be in-depth analyses

and comparison of differences between ground-based estimates of soil moisture, MOD16 data from this study, original MOD16 AET and PET, and ground-observed vegetation characteristics. Similarly, a better spatially varying estimated value for albedo could improve ET estimation. Satellite data are considered to be of lower quality than groundbased estimates (see Section 1). The ratio of uncertainty of satellite data to ground-based estimates was therefore chosen as a conservative 1.5. This estimate is debatable, as the exact figure is not known. The reason for that is that the uncertainty in global meteorology data, cloud cover, land cover, vegetation characteristics are unknown. Further research of the combined effect of cloud cover and vegetation characteristics, and the difference of global with national meteorological data could help to better quantify this, and is a recommendation for further research. 6. Conclusions As global satellite data have potential for use on regional or catchment scale, MODIS MOD16 data have been tested to improve ET estimation on a 1 km × 1 km monthly resolution in New Zealand. While existing techniques consider ground-based estimates as the ‘gold standard’, this research incorporates uncertainty of both groundobserved and satellite ET. Uncertainty analysis performed on groundobserved ET0 shows that its uncertainty varies between 10% and 40% of ET0, for high and low ET0, respectively. It also shows that ET0 is most sensitive to temperature, followed by solar radiation, cloudiness ratio and relative humidity. Using this analysis, a set of correlated random variables, and a Monte-Carlo fitting approach, satellite data becomes a ‘soft interpolator’ between the ground-based estimates. The resulting ET0 estimates also contain an uncertainty estimate. The proposed method enhances the capability of using global satellite data products on a catchment scale, hereby abiding the local used Penman standard.

110

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

Fig. 9. Monthly MOD160 (red) and ground-based ET0 estimates (in black) in January; April; July; and October of 2010.

Comparison between the P–M derived MOD16 PET and groundobserved ET0 suggests that it is very well possible that the original MOD16 data can be directly applied in a PET and AET estimation. However, quantitative data to further validate this are lacking. Recommendations for further research are therefore further quantitative analyses. Also, temporal up-scaling to daily estimates with existing gridded (VCS) data is recommended. Acknowledgements This research has been performed in the Smart Aquifer Characterisation Programme, funded by the Ministry of Business, Innovation and Employment, New Zealand. It is part of a PhD study of the author at the University of Waikato, New Zealand. This project has received co funding from the European Union's Seventh Programme for research technological development and demonstration under grant agreement No 603608, eartH2Observe. I would like to thank Paul White, Mike Friedel (GNS Science), and Jaap Schellekens (Deltares) for inspiring me with ideas on error propagation and terminology of ET. I thank the reviewers for their in-depth review, which has helped improve my research. I would furthermore like to thank Andrew Tait and MS Srinivasan from NIWA, New Zealand, for generously sharing ground-observed and gridded VCS ET and its background information. Last but not least, I would like to thank Qiaozhen Mu, University

of Montana, for supplying me with MOD16 data, their updates, and further information. Appendix A. Calculation of input component in Penman and Penman–Monteith equations The value of es [kPa] can be approximated as a function of air temperature only through (Lin & Hubbard, 2004; World Meteorological Organization, 2008):  es ¼ 0:6108 exp

 17:27T : 237:3 þ T

ðA:1Þ

If es is known, ea can be derived by measuring RH: RH ¼ 100

ea : es

ðA:2Þ

The gradient of saturation vapour pressure over temperature Δ is a function of es and T (Murray, 1967; Tetens, 1930): Δ¼

des 4098es ¼ : dT ð237:3 þ T Þ2

ðA:3Þ

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

Fig. 10. Monthly ET0 for all months at 6 randomly picked climate stations for MOD160 (red dots) and ground-based ET0 (black triangles).

Fig. 11. Mean annual MOD160 ET0 (2000–2014).

Fig. 12. Mean annual σET0 of 2000–2014 MOD160 ET0.

111

112

R.S. Westerhoff / Remote Sensing of Environment 169 (2015) 102–112

Rn is the difference between the net incoming short-wave radiation from the sun at the earth's surface Sn and outgoing long-wave radiation (from the earth) Ln: Rn ¼ Sn −Ln :

ðA:4Þ

Both Ln and Sn can be estimated by empirical relations (Hendriks, 2010): pffiffiffiffiffi  n Ln ¼ 4:903  10−9 ðT þ 273:2Þ4 ð0:34−0:14 ea Þ 0:25 þ 0:75 ; ðA:5Þ N  n Sn ¼ ð1−α Þ 0:25 þ 0:50 S0 ; N

ðA:6Þ

where: n/N is the ratio of bright sunshine hours per day n and the total possible hours of sunshine per day N; α is the earth's albedo (solar reflection coefficient); and S0 is the short-wave radiation at the top of the atmosphere. Sn relates to measured solar radiation Rad at the surface as (Burman & Pochop, 1994; NIWA, 2014a, 2014b). Sn ¼ ð1−α ÞRad

ðA:7Þ

References Ahrens, C.D. (2007). Meteorology Today: An Introduction to Weather, Climate, and the Environment (8th ed.). Belmont, CA: Thomson/Brooks/Cole. Allen, R.G. (1998). Crop evapotranspiration: Guidelines for computing crop water requirements. Number 56 in FAO Irrigation and Drainage Paper. Rome: Food and Agriculture Organization of the United Nations. 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. Journal of Geophysical Research, 112. http://dx.doi.org/10.1029/2006JD007506. Australian National University (2014). Anusplin vrsn 4.4. http://fennerschool.anu.edu.au/ research/products/anusplin-vrsn-44 (accessed October 31, 2014) Bastiaanssen, W., Pelgrum, H., Wang, J., Ma, Y., Moreno, J., Roerink, G., et al. (1998). A remote sensing surface energy balance algorithm for land (SEBAL). part 2: Validation. Journal of Hydrology, 212-213, 213–229. http://dx.doi.org/10.1016/S00221694(98)00254-6. Burman, R.D., & Pochop, L.O. (1994). Evaporation, evapotranspiration and climatic data. Number 22 in Developments in Atmospheric Science. New York: Elsevier, Amsterdam. Chirlin, G.R., & Wood, E.F. (1982). On the relationship between kriging and state estimation. Water Resources Research, 18, 432–438. http://dx.doi.org/10.1029/ WR018i002p00432. Coleman, G., & DeCoursey, D.G. (1976). Sensitivity and model variance analysis applied to some evaporation and evapotranspiration models. Water Resources Research, 12, 873–879. eartH2Observe (2015). earth2observe website. http://www.earth2observe.eu (accessed February 18, 2015) Estévez, J., Gavilàn, P., & Berengena, J. (2009). Sensitivity analysis of a Penman–Monteith type equation to estimate reference evapotranspiration in southern Spain. Hydrological Processes, 23, 3342–3353. GEO (2014a). The GEOSS Water Strategy: From Observations to Decisions. Tokyo: Japan Aerospace Exploration Agency, Genveva (URL: ftp://ftp.earthobservations.org/ TEMP/Water/GEOSS_WSR_Full_Report.pdf. accessed December 12, 2014). GEO (2014b). Group on earth observation. http://www.earthobservations.org (accessed December 10, 2014) GEWEX (2014). Global energy and water exchange project. http://www.gewex.org (accessed December 10, 2014) GLOWASIS (2015). GLOWASIS website. http://glowasis.eu (accessed February 18, 2015) Gong, L., Xu, C.Y., Chen, D., Halldin, S., & Chen, Y.D. (2006). Sensitivity of the Penman– Monteith reference evapotranspiration to key climatic variables in the Changjiang (Yangtze river) basin. Journal of Hydrology, 329, 620–629. Hendriks, M.R. (2010). Introduction to Physical Hydrology. New York: Oxford University Press, Oxford. Howard, K., & Lloyd, J. (1979). The sensitivity of parameters in the Penman evaporation equations and direct recharge balance. Journal of Hydrology, 41, 329–344. Irmak, A. (Ed.). (2012). Evapotranspiration — Remote Sensing and Modeling (Hardcover ed.). InTech. http://dx.doi.org/10.5772/725. Leathwick, J.R., & Stephens, R.T.T. (1998). Climate surfaces for New Zealand. Landcare Research Contract Report LC9798, 126. (pp. 119). Li, Z., Cribb, M., & Chang, F.L. (2005). Natural variability and sampling errors in solar radiation measurements for model validation over the atmospheric radiation measurement southern Great Plains region. Journal of Geophysical Research, 110. http://dx. doi.org/10.1029/2004JD005028. Lin, X., & Hubbard, K.G. (2004). Uncertainties of derived dewpoint temperature and relative humidity. Journal of Applied Meteorology, 43, 821–825. http://dx.doi.org/10.1175/ 2100.1.

Mathworks (2014a). Multivariate normal distribution – MATLAB & simulink – MathWorks Australia. http://www.mathworks.com.au/help/stats/multivariatenormal-distribution.html (accessed September 22, 2014) Mathworks (2014b). Smoothing splines – MATLAB & Simulink – MathWorks Australia. http://www.mathworks.com.au/help/curvefit/smoothing-splines.html (accessed September 1, 2014) McMahon, T.A., Peel, M.C., Lowe, L., Srikanthan, R., & McVicar, T.R. (2013). Estimating actual, potential, reference crop and pan evaporation using standard meteorological data: A pragmatic synthesis. Hydrology and Earth System Sciences, 17, 1331–1363. http://dx.doi.org/10.5194/hess-17-1331-2013. Miralles, D., De Jeu, R., Gash, J., Holmes, T., & Dolman, A. (2011). Magnitude and variability of land evaporation and its components at the global scale. Hydrology and Earth System Sciences, 15, 967–981. http://dx.doi.org/10.5194/hess-15-967-2011. Miralles, D.G., Holmes, T.R.H., De Jeu, R.A.M., Gash, J.H., Meesters, A.G.C.A., & Dolman, A.J. (2011). Global land-surface evaporation estimated from satellite-based observations. Hydrology and Earth System Sciences, 15, 453–469. http://dx.doi.org/10.5194/hess-15453-2011. Monteith, J. (1965). Evaporation and environment. Symposium of the Society of Experimental Biology, 19, 205–224. Mu, Q., Heinsch, F., Zhao, M., & Running, S. (2007). Development of a global evapotranspiration algorithm based on MODIS and global meteorology data. Remote Sensing of Environment, 111, 519–536. Mu, Q., Zhao, M., & Running, S. (2011). Running, improvements to a MODIS global terrestrial evapotranspiration algorithm. Remote Sensing of Environment, 115, 1781–1800. http:// dx.doi.org/10.1016/j.rse.2011.02.019. Murray, F. (1967). On the computation of saturation vapor pressure. Journal of Applied Meteorology, 6, 203–204. NIWA (2014a). Cliflo. evaporation: Methodology. http://cliflo-niwa.niwa.co.nz/pls/niwp/ wh.do_help?id=ls_evap_method (accessed November 11, 2014) NIWA (2014b). CliFlo: the national climate database. http://cliflo.niwa.co.nz/ (accessed January 12, 2015) NTSG (2013). MODIS global evapotranspiration project (MOD16). http://www.ntsg.umt. edu/project/mod16 (accessed January 14, 2014) Oke, T.R. (1992). Boundary Layer Climates. New York: Routledge, London. Penman, H. (1948). Natural evaporation from open water, bare soil and grass. Proceedings of the Royal Society, Series A, 193, 120–145. Penman, H. (1963). Vegetation and hydrology. Tech. Comm., vol. 53, Harpenden, England: Commonwealth Bureau of Soils. Rose, C., & Smith, M. (1996). The multivariate normal distribution. The Mathematica Journal, 6, 32–37. Saxton, K.E. (1975). Sensitivity analyses of the combination evapotranspiration equation. Agricultural Meteorology, 15, 343–353. http://dx.doi.org/10.1016/00021571(75)90031-X. Su, Z. (2002). The surface energy balance system (SEBS) for estimation of turbulent heat fluxes. Hydrology and Earth System Sciences, 6, 85–100. http://dx.doi.org/10.5194/ hess-6-85-2002. Tait, A. (2014). Personal Communication. Principal Scientist at the National Climate Centre, NIWA. Tait, A., & Woods, R. (2007). Spatial interpolation of daily potential evapotranspiration for New Zealand using a spline model. Journal of Hydrometeorology, 8, 430–438. http:// dx.doi.org/10.1175/JHM572.1. Tait, A., & Zheng, X. (2007). Analysis of the spatial interpolation error associated with maps of median annual climate variables. Technical Report. NIWA (URL: https:// www.niwa.co.nz/sites/niwa.co.nz/files/import/attachments/Climate_Maps_Error_ Analysis.pdf. 21 pages, accessed July 7, 2014). Tait, A., Henderson, R., Turner, R., & Zheng, X. (2006). Thin plate smoothing spline interpolation of daily rainfall for New Zealand using a climatological rainfall surface. International Journal of Climatology, 26, 2097–2115. http://dx.doi.org/10.1002/joc.1350. Taylor, J.R. (1997). An Introduction to Error Analysis: The Study of Uncertainties in Physical Measurements (2nd ed.). Sausalito, Calif: University Science Books. Tellinghuisen, J. (2001). Statistical error propagation. The Journal of Physical Chemistry. A, 105, 3917–3921. http://dx.doi.org/10.1021/jp003484u. Tetens, O. (1930). Über einige meteorologische begriffe. Zeitschrift für Geophysik, 6, 297–309. Thornthwaite, C. (1948). An approach toward a rational classification of climate. Geographical Review, 38, 55–94. Thornthwaite, C., & Mather, J. (1957). Instructions and tables for computing potential evapotranspiration and the water balance, Centerton, N.J., Laboratory of Climatology. Publications in Climatology, 10, 185–311. UCAR (2014). Unidata | NetCDF. http://www.unidata.ucar.edu/software/netcdf/ (accessed 1 July, 2014) Van der Tol, C., & Parodi, G. (2012). Guidelines for remote sensing of evapotranspiration. Evapotranspiration — Remote Sensing and Modeling. InTech, 227–250. Verhoef, A., & Egea, G. (2014). Modeling plant transpiration under limited soil water: Comparison of different plant and soil hydraulic parameterizations and preliminary implications for their use in land surface models. Agricultural and Forest Meteorology, 191, 22–32. Wilson, K., Goldstein, A., Falge, E., Aubinet, M., Baldocchi, D., Berbigier, P., et al. (2002). Energy balance closure at FLUXNET sites. Agricultural and Forest Meteorology, 113, 223–243. Woods, R., Hendrikx, J., Henderson, R., & Tait, A. (2006). Estimating mean flow of New Zealand rivers. Journal of Hydrology. New Zealand, 45, 95–110. World Meteorological Organization (2008). Guide to Meteorological Instruments and Methods of Observation. Geneva, Switzerland: World Meteorological Organization. Zhao, M., Running, S.W., & Nemani, R.R. (2006). Sensitivity of Moderate Resolution Imaging Spectroradiometer (MODIS) terrestrial primary production to the accuracy of meteorological reanalyses. Journal of Geophysical Research, 111. http://dx.doi.org/10. 1029/2004JG000004.