Assessing Crop Coefficients for Natural Vegetated Areas Using Satellite Data and Eddy Covariance Stations

The Food and Agricultural Organization (FAO) method for potential evapotranspiration assessment is based on the crop coefficient, which allows one to relate the reference evapotranspiration of well irrigated grass to the potential evapotranspiration of specific crops. The method was originally developed for cultivated species based on lysimeter measurements of potential evapotranspiration. Not many applications to natural vegetated areas exist due to the lack of available data for these species. In this paper we investigate the potential of using evapotranspiration measurements acquired by micrometeorological stations for the definition of crop coefficient functions of natural vegetated areas and extrapolation to ungauged sites through remotely sensed data. Pastures, deciduous and evergreen forests have been considered and lower crop coefficient values are found with respect to FAO data.


Introduction
Accurate assessment of evapotranspiration is fundamental for efficient water resource management, for the design and management of water supply reservoirs, for designing and scheduling irrigation systems, and for environmental assessment [1]. Moreover, trends of modifications of water availability around the world, have prompted scientific investigations of the effect of climate change on the hydrological cycle [2][3][4][5]. In order to study the effects of climate change for several years in advance, the use of parsimonious models for the assessment of evapotranspiration is needed in order to maintain the simulation time within acceptable limits and because of limited availability of meteorological data at a spatial and temporal resolution sufficient to accurately capture the dynamics of hydrological processes.
Latent heat flux (LE), i.e., evapotranspiration (ET), is the joint term that rules the relationship between water and energy fluxes through vegetated soil and shallow atmosphere layer, since it is the only term present both in the water as well as in the energy balance. Despite its important role, direct measurements of ET are neither frequent nor simple, due to the complexity of the process, its spatial scale of representativeness and plant phenology. These issues often affect the reliability of direct measures of water losses, so that the scientific community proposed in the late forties onwards several methods to estimate ET based on meteorological measures [6]. A simplified method is based on the crop coefficient (kc) which embodies all the physiologic characteristics of a specific plant allowing one to relate a reference evapotranspiration of well irrigated grass to the potential evapotranspiration of (a) The definition of crop coefficient curves for natural area derived from eddy covariance data to be used in hydrological modelling to compute effective evapotranspiration. (b) Assessing the reliability and potentiality of using satellite data and conventional meteorological measurements for crop coefficient estimates in natural vegetated areas where eddy covariance stations are not available.
In this paper the proposed procedure is applied for three types of natural coverage: pastures, deciduous and evergreen forests.

Materials and Methods
Much literature discusses "potential" and "actual" ET, the former being the maximum evapotranspiration which would occur if no factors other than energy supply and atmospheric demand limiting the ET rate [32]. Conversely, "actual" ET is the amount of water evapotranspirated in the case that water availability becomes a limiting factor [33].
Potential evapotranspiration can be derived from ground measurements, lysimeter or eddy covariance stations, or computed using a wide range of equations that are very well assessed in the research community, such as the Prietsley-Taylor equation [34], Hargreaves [8] or Penman Monteith [6]. In this study the Penman-Monteith equation has then been applied.

Penman-Monteith Equation
Penman-Monteith's equation, that combines the energy balance with the mass transfer method, allows computing a potential evapotranspiration (ETP PenMon (mm·day −1 ) as: where R n is the net radiation (MJ·m −2 ·d −1 ), G is the soil round heat flux (MJ·m −2 ·d −1 ), (e 0 − e a ) is the vapour pressure deficit of the air (kPa), ρ a is the mean air density (kg·m −3 ), c p is the specific heat of the air (MJ·kg −1 · • C −1 ), ∆ is the slope of the saturation vapour pressure (kPa· • C −1 ), γ is the psychrometric constant (kPa· • C −1 ) and r c is the crop resistance (s·m −1 ) and r a is the aerodynamic resistance (s·m −1 ). In particular net radiation is computed as sum of longwave and shortwave radiation: where α is albedo, R s is the incoming shortwave radiation (MJ·m −2 ·d −1 ) and R nl is the net longwave radiation (MJ·m −2 ·d −1 ). The aerodynamic resistance is a function of wind velocity (u) (m·s −1 ) and canopy height (h c ) (m) [35]: where z u and z rh are respectively the measurement heights of wind velocity and relative air humidity (m). The crop resistance is computed as [36]: where r smin is the minimum stomatal resistance (s·m −1 ) and LAI is leaf area index. G is approximated as a fraction (0.1 during daytime [7]) of net radiation.

Sensitivity Analysis
A sensitivity analysis on kc is performed to understand how much it deviates from its mean value. A sensitivity coefficient is computed in respect to kc for each parameter included in crop coefficient estimates (e.g., LAI, albedo and vegetation height) and then is adimensionalized so that the results are comparable. The coefficient of sensitivity is computed as: where S Kc is the sensitivity coefficient and V i is the ith parameter. A positive value of S Kc means that the parameter will increase as kc increases and the opposite behaviour for negative values of S Kc .

Eddy Covariance Technique
The eddy covariance method determines the surface fluxes as the sum of turbulent fluxes, measured above the surface, and the flux divergence between the surface and the eddy covariance measurement level [37]. The basic equations to estimate latent and sensible heat fluxes are: where λ is the vaporization latent heat, ρ the air density and w H 2 O the covariance between vertical wind velocity and concentration of vapor in the air. C p is the specific heat at constant pressure and w T is the covariance between vertical wind velocity and temperature. The symbol ( ) indicates instantaneous fluctuation from the time-averaged values of a specific variable in according with Reynolds decomposition of a meteorological stochastic signal [38]. The reliability of flux measurements depends on different theoretical assumptions of the eddy covariance technique [17,39], among which the most important are the horizontal homogeneity, the stationarity and the mean vertical wind speed equal to zero during the averaging period. A number of corrections needs to be applied to obtain high quality fluxes flowing procedures which are now well assessed in literature [40,41].

The FAO Crop Coefficient
According to FAO Paper no. 56 [7], the crop coefficient is given by: where ET 0 is the reference potential evapotranspiration (mm·day −1 ) and ETP is a generic potential evapotranspiration (mm·day −1 ). According to FAO Paper no. 56, three values of kc are defined over time in accordance to the crop development stage: (i) initial stage (kc -ini ), (ii) mid-season stage (kc -mid ) and the harvest stage (kc -end ). The reference evapotranspiration, ET 0 , is a potential evapotranspiration which derives from Penman-Monteith's equation and is defined for a reference surface defined as an "hypothetical crop with an assumed height of 0.12 m having a surface resistance of 70 s·m −1 and an albedo of 0.23, closely resembling the evaporation of an extension surface of green grass of uniform height, actively growing and adequately watered" [7].
ET 0 (mm·day −1 ) is computed as: where u 2 is wind speed at 2 m above ground surface (m·s −1 ) expressed as: For natural vegetation [7] proposed to assess crop coefficient as a function of leaf area index (LAI) where: where u 2 is the mean value for wind speed at 2 m height during the mid-season (m·s −1 ), RH min is the mean value for minimum daily relative humidity during mid-season (%), and h is the mean maximum plant height (m).

Crop Coefficients
In this paper, crop coefficient is assessed in three ways that are different as regards ETP estimates: kc PenMon,sat = ETP PenMon,sat ET 0 (15) kc PenMon,ground = ETP PenMon,ground ET 0 (16) where ETP eddy is the potential evapotranspiration derived by eddy covariance stations, ETP PenMon,sat is the potential evapotranspiration estimated with Penman-Monteith's equation using leaf area index and albedo retrieved from remotely sensed images and ground measured meteorological forcings, and ETP PenMon,ground is the potential evapotranspiration estimated with Penman-Monteith's equation using data coming from ground measurements.
Since the crop coefficient approach proposed by FAO is intended for potential evapotranspiration estimate, kc eddy is significant only when vegetation can evapotranspirate with no limitation in unlimited water conditions. As it can be expected that, when available moisture of the soil is enough, ET from eddy covariance station should be similar to the one calculated with the Penman-Monteith's equation (Equation (1)), in the present work kc eddy is computed only when the condition ETP eddy ETP PenMon,ground ≈ 1 (±0.05) is satisfied.
In addition to above mentioned methods for crop coefficient assessment, Equation (12) is applied (kc FAO ). A schematic of the methodology followed in this research is illustrated in Figure 1.   (16) where ETPeddy is the potential evapotranspiration derived by eddy covariance stations, ETPPenMon,sat is the potential evapotranspiration estimated with Penman-Monteith's equation using leaf area index and albedo retrieved from remotely sensed images and ground measured meteorological forcings, and ETPPenMon,ground is the potential evapotranspiration estimated with Penman-Monteith's equation using data coming from ground measurements.
Since the crop coefficient approach proposed by FAO is intended for potential evapotranspiration estimate, kceddy is significant only when vegetation can evapotranspirate with no limitation in unlimited water conditions. As it can be expected that, when available moisture of the soil is enough, ET from eddy covariance station should be similar to the one calculated with the Penman-Monteith's equation (Equation (1) In addition to above mentioned methods for crop coefficient assessment, Equation (12) is applied (kcFAO). A schematic of the methodology followed in this research is illustrated in Figure 1.

Sites and Data
Monitoring took place in four sites representative of three natural ecosystems: the first one is located in a grass field in the mountain area of Torgnon in Italy [45 •  The study site of Torgnon [42,43] is an abandoned pasture located in northwestern Italian Alps (Aosta Valley, Italy) at an elevation of 2160 m a.s.l. The site is part of the FLUXNET network. This site is characterized by an intra-Alpine semi-continental climate, with an annual mean temperature of 3.1 • C and a mean annual precipitation of about 880 mm. From the end of October to late May, the site is covered by a thick snowcover (90-120 cm) which limits the growing period, and hence plants evapotranspiration, to an average of five months. Dominant vegetation is composed by matgrass (Nardus stricta L., Festuca nigrescens All., Arnica montana L., Carex sempervirens Vill., Geum montanum L., Anthoxanthum alpinum L.L., Potentilla aurea L., Trifolium alpinum L..) and the soil is classified as Cambisol (WRB classification).
The eddy covariance method is used to measure the fluxes of H 2 O between the ecosystem and the atmosphere. Measurement of wind speed in the three components is performed by a CSAT3 three-dimensional sonic anemometer (Campbell Scientific, Logan, UT, USA), while H 2 O vapor air densities were measured by a LI-7500 open-path infrared gas analyzer (LICOR, Lincoln, NE, USA). Instruments were placed at 2.5 m above the ground. Data from both the anemometer and the IRGA were measured at 10 Hz. A weather station is installed close to the eddy covariance tower to continuously measure additional meteorological variables. Air and soil temperatures are measured respectively by a HMP45 (Vaisala, Helsinki, Finland) and with temperature probes type therm107 (Campbell Scientific) at different depths. Soil water content is assessed with soil water reflectometers model CS616 (Campbell Scientific). The variables of interest for this study, i.e., soil heat flux and net radiation are measured respectively by HFP01 plates (Hukseflux, Delft, The Netherlands) and by a CNR4 (Kipp & Zonen, Delft, The Netherlands) net radiometer. Snow height is measured with a SR50A sonic snow depth sensor (Campbell Scientific).
Since the grassland is unmanaged, canopy structure varies during the season following the phenological development. Measurements of canopy height and Leaf Area Index (LAI) at site are available for the study years, in detail: canopy height is measured and phytomass is sampled every two weeks at 12 selected plots of 40 cm × 40 cm in the study area. Total LAI and canopy height for the 12 plots are averaged to obtain a site-level value.

Satellite Data
Leaf area index and albedo are retrieved from satellite data measured from MODIS radiometers on board TERRA and AQUA NASA satellites (https://earthdata.nasa.gov/). These satellites fly over the same area twice a day, providing a night and a day image. The values of LAI [48] and albedo [49] are available products which are averaged over a range of eight days to avoid problems of cloudy conditions. These parameters can be considered constant during this period of time.

Results and Discussion
kc values are computed at daily scale from measured eddy covariance data at 30 min or 1 h. A day of data is considered only if at least the 80% of daytime data is available (when R n > 0), otherwise the entire day is deleted. Other gaps are also present in the observed series due to the identification of potential evapotranspiration conditions. In fact days characterized by soil water stress conditions are eliminated from the data series.
Finally in order to compute the difference between the computed crop coefficient values, RMSE is evaluated as: where X sim ith is the ith simulated variable, X obs ith is the ith measured variable, n the sample size, and X obs the average observed variable.

Results of the Sensitivity Analysis
In Figure 2 the sensitivity coefficients for kc are reported for LAI, albedo and vegetation height considering a variation from the mean value of ±50%, using the meteorological variables to compute energy fluxes measured by an eddy covariance station located in a maize field in Northern Italy [50]. +50%. Instead the crop coefficient is highly influenced by the variability of leaf area index and in particular decreases in exponential way from −50% variation of LAI with Skc equal to 0.44 till values of 0.13 with a variation of +50% of LAI, so that the crop coefficient will be more sensitive during the growing period of the vegetation. These results are also confirmed by [51] who found that the most sensitive parameter in ET0 estimation is rc.

Pasture: Torgnon
In Figure Table 1, while in Table 2 the lengths of the growth stages are reported.
During the year 2009, two crop coefficient values from eddy covariance measurements can be identified: the first one is relative to a completely developed grass with a mean values of 0.8, while the second one is relative to the senescence phase when grass starts to become yellow with a mean value of 0.47. A similar distribution is found also during 2010, with the first kceddy equal to 0.88 and the second to 0.43. Albedo and vegetation height have a low influence on crop coefficient estimates; in fact S Kc remains almost constant and equal to 0.1 for variations of albedo and vegetation height from −50% to +50%. Instead the crop coefficient is highly influenced by the variability of leaf area index and in particular decreases in exponential way from −50% variation of LAI with S Kc equal to 0.44 till values of 0.13 with a variation of +50% of LAI, so that the crop coefficient will be more sensitive during the growing period of the vegetation. These results are also confirmed by [51] who found that the most sensitive parameter in ET 0 estimation is r c .

Pasture: Torgnon
In Figure Table 1, while in Table 2 the lengths of the growth stages are reported.
During the year 2009, two crop coefficient values from eddy covariance measurements can be identified: the first one is relative to a completely developed grass with a mean values of 0.8, while the second one is relative to the senescence phase when grass starts to become yellow with a mean value of 0.47. A similar distribution is found also during 2010, with the first kc eddy equal to 0.88 and the second to 0.43.
When the crop coefficient is computed with Penman-Monteith's equation, a kc PenMon,ground of 0.76 is obtained during 2010 for the first period and 0.39 for the senescence period. These values are quite similar to those obtained from eddy covariance data, showing that for pasture correct kc values can be defined using meteorological data from basic network stations which are more available than from eddy covariance towers. +50%. Instead the crop coefficient is highly influenced by the variability of leaf area index and in particular decreases in exponential way from −50% variation of LAI with Skc equal to 0.44 till values of 0.13 with a variation of +50% of LAI, so that the crop coefficient will be more sensitive during the growing period of the vegetation. These results are also confirmed by [51] who found that the most sensitive parameter in ET0 estimation is rc.

Pasture: Torgnon
In Figure 3 the computed crop coefficient values are reported for 2009 and 2010. The selected data belong to snow-free periods from day 183 to 305 for 2009 and from 142 to 302 for 2010. Results for kc are summarized in Table 1, while in Table 2 the lengths of the growth stages are reported.
During the year 2009, two crop coefficient values from eddy covariance measurements can be identified: the first one is relative to a completely developed grass with a mean values of 0.8, while the second one is relative to the senescence phase when grass starts to become yellow with a mean value of 0.47. A similar distribution is found also during 2010, with the first kceddy equal to 0.88 and the second to 0.43. The analyses are then performed using satellite information of LAI and kc PenMon,sat are equal to 0.76 for the first period and 0.49 for the second period during 2009, while 0.76 and 0.36 for 2010. kc PenMon,sat are similar to kc PenMon,ground and kc eddy , showing that this technique can be applied in a distribute way in river basins coupling basic network stations for meteorological data and satellite data for vegetation information.
Finally, RMSE values have been computed between kc PenMon,sat and kc eddy and RMSE values equal to 0.15 and 0.19 are found for 2009 and 2010 respectively. When kc PenMon,ground and kc eddy are compared for 2010, a RMSE of 0.18 is obtained.

Chestnut Ridge
In Figure 4 crop coefficient values are reported for the deciduous forest of Chestnut Ridge. The data belong to the whole year, so that the complete seasonal cycle can be analyzed differencing the behaviors during periods with or without leaves (Table 1). In Table 2 the lengths of the growth stages are reported.
Sensors 2017, 17, 2664 9 of 14 be defined using meteorological data from basic network stations which are more available than from eddy covariance towers. The analyses are then performed using satellite information of LAI and kcPenMon,sat are equal to 0.76 for the first period and 0.49 for the second period during 2009, while 0.76 and 0.36 for 2010. kcPenMon,sat are similar to kcPenMon,ground and kceddy, showing that this technique can be applied in a distribute way in river basins coupling basic network stations for meteorological data and satellite data for vegetation information.
Finally, RMSE values have been computed between kcPenMon,sat and kceddy and RMSE values equal to 0.15 and 0.19 are found for 2009 and 2010 respectively. When kcPenMon,ground and kceddy are compared for 2010, a RMSE of 0.18 is obtained.

Chestnut Ridge
In Figure 4 crop coefficient values are reported for the deciduous forest of Chestnut Ridge. The data belong to the whole year, so that the complete seasonal cycle can be analyzed differencing the behaviors during periods with or without leaves (Table 1). In Table 2 the lengths of the growth stages are reported. During the winter period, kceddy has low values equal to 0.19. Similar results are obtained during the autumn season when the leaves are not present as well. kc from eddy covariance tower is then equal to 0.2. When these periods are analyzed using satellite information, much lower values are obtained with kcPenMon,sat for the initial period equal to 0.04 and kc-end to 0.05 for both years. This behaviour is linked to values near zero or zero of LAI from satellite during these winter-autumn periods.
The growth and decrease branches of the two kc curves are quite in accordance, according to Figure 4, as well as for the summer period when kc-mid from eddy covariance station is equal to 0.47 and from satellite data to 0.48. Overall, analyzing the estimates accuracy for the whole year, RMSE obtained from the comparison between kc computed from eddy covariance data and from satellite data is found to be equal to 0.18.

Duke Forest
Crop coefficient results for the years 2001 and 2002 are reported in Figure 5. kceddy during 2001 is low during winter with a mean value of 0.11 while kcPenMon,sat is 0.07. Similar values are obtained for 2002 with kcPenMon,sat being 0.04 and kceddy 0.09. These results are obtained also during autumn (Table 1). During the summer season, when trees are completely lush, as expected, crop coefficient values are During the winter period, kc eddy has low values equal to 0.19. Similar results are obtained during the autumn season when the leaves are not present as well. kc from eddy covariance tower is then equal to 0.2. When these periods are analyzed using satellite information, much lower values are obtained with kc PenMon,sat for the initial period equal to 0.04 and kc -end to 0.05 for both years. This behaviour is linked to values near zero or zero of LAI from satellite during these winter-autumn periods.
The growth and decrease branches of the two kc curves are quite in accordance, according to Figure 4, as well as for the summer period when kc -mid from eddy covariance station is equal to 0.47 and from satellite data to 0.48. Overall, analyzing the estimates accuracy for the whole year, RMSE obtained from the comparison between kc computed from eddy covariance data and from satellite data is found to be equal to 0.18.  (Table 2).

Intercomparison
A comparison between Chestnut Ridge and Duke Forests results is performed in order to define a general crop coefficient curve which is representative for deciduous forests. According to Table 1, results of crop coefficient from eddy covariance data are quite similar between the two stations for the initial stage value for all the analysed year ranging between 0.09 and 0.19, as well as for the final stage with values between 0.12 and 0.2. kc-mid is also quite defined with values between 0.43 and 0.48.

Evergreen Forest: Black Hills
In Figure 6 crop coefficient values are reported for 2006 and 2007 for the Black Hills station. Even though the station is located in an evergreen forest, evapotranspiration substantially decreases during the winter period and so the crop coefficient does. In Table 1 kceddy has low values equal to 0.05 during the winter period for both years. Similar valus are obtained for kcPenMon,sat. The final stage has similar low values equal to 0.03 and 0.07 during 2006 from eddy covariance data and from satellite data respectively, while for 2007 kceddy is equal to 0.03 while kcPenMon,sat to 0.04.
The growth and decrease branches of the two kc curves are quite in accordance, according to Figure 4 and Table 2, as well as for the summer period when kc-mid from eddy covariance station is

Intercomparison
A comparison between Chestnut Ridge and Duke Forests results is performed in order to define a general crop coefficient curve which is representative for deciduous forests. According to Table 1, results of crop coefficient from eddy covariance data are quite similar between the two stations for the initial stage value for all the analysed year ranging between 0.09 and 0.19, as well as for the final stage with values between 0.12 and 0.2. kc -mid is also quite defined with values between 0.43 and 0.48.

Evergreen Forest: Black Hills
In Figure 6 crop coefficient values are reported for 2006 and 2007 for the Black Hills station. Even though the station is located in an evergreen forest, evapotranspiration substantially decreases during the winter period and so the crop coefficient does.

Intercomparison
A comparison between Chestnut Ridge and Duke Forests results is performed in order to define a general crop coefficient curve which is representative for deciduous forests. According to Table 1, results of crop coefficient from eddy covariance data are quite similar between the two stations for the initial stage value for all the analysed year ranging between 0.09 and 0.19, as well as for the final stage with values between 0.12 and 0.2. kc-mid is also quite defined with values between 0.43 and 0.48.

Evergreen Forest: Black Hills
In Figure 6 crop coefficient values are reported for 2006 and 2007 for the Black Hills station. Even though the station is located in an evergreen forest, evapotranspiration substantially decreases during the winter period and so the crop coefficient does. In Table 1 kceddy has low values equal to 0.05 during the winter period for both years. Similar valus are obtained for kcPenMon,sat. The final stage has similar low values equal to 0.03 and 0.07 during 2006 from eddy covariance data and from satellite data respectively, while for 2007 kceddy is equal to 0.03 while kcPenMon,sat to 0.04.
The growth and decrease branches of the two kc curves are quite in accordance, according to Figure 4 and Table 2, as well as for the summer period when kc-mid from eddy covariance station is In Table 1 kc eddy has low values equal to 0.05 during the winter period for both years. Similar valus are obtained for kc PenMon,sat . The final stage has similar low values equal to 0.03 and 0.07 during 2006 from eddy covariance data and from satellite data respectively, while for 2007 kc eddy is equal to 0.03 while kc PenMon,sat to 0.04.
The growth and decrease branches of the two kc curves are quite in accordance, according to Figure 4 and Table 2, as well as for the summer period when kc -mid from eddy covariance station is equal to 0.17 and from satellite data to 0.2 for 2006 and to 0.18 and 0.2 respectively for 2007. Overall, analyzing the estimates accuracy for the whole year, RMSE between kc from eddy covariance data and from satellite data is equal to 0.07 in 2006 and to 0.05 in 2007.

Conclusions
In this paper, the potentiality of using evapotranspiration measurements from micro-meteorological stations for definition of crop coefficient functions of natural vegetated area has been evaluated. Sensitivity analysis shows that albedo and vegetation height have a low influence on crop coefficient estimates, while these are highly influenced by the variability of leaf area index.
For pasture, crop coefficient obtained from eddy covariance data are similar to the values obtained with ground measured data and those obtained using leaf area index and albedo retrieved from remotely sensed images. This leads to two important results: (1) measurements from widespread available stations provide robust estimate of crop coefficient, and (2) satellite data are powerful to compute crop coefficient in a distribute way in river basins wherever basic ground measured meteorological data are available.
The potentials of remotely sensed data have been confirmed using data available from the FLUXNET network both on deciduous and evergreen forests. The crop coefficient values for the considered coverage of pastures, deciduous and evergreen forests are found to be lower than the FAO values during the mid-season stage, showing the necessity to check FAO data before applying them to areas with climatic and vegetation characteristics different form the ones considered by FAO.