Estimation of ecosystem evapotranspiration in a Robinia pseudoacacia L. plantation with the use of the eddy covariance technique and modeling approaches Open Access

Evapotranspiration (ET) is one of the greatest components of energy fluxes and water cycling in the atmosphere (Katul et al. 2012). It is responsible for the returning of up to 2/3 of precipitation water back to the atmosphere on average globally (Jasechko et al. 2013; Wei et al. 2017; Yao et al. 2017), an amount that can be as high as 95% in dry ecosystems (Raz-Yaseef et al. 2010; Stoy et al. 2019). The accurate estimation of evapotranspiration is essential for a better understanding of the global water cycle (Cheng et al. 2011) and the improvement of water resource management (Mu et al. 2013a; Paschalis et al. 2018).

For ET estimation, several methods have been developed, including the use of the eddy covariance method (Baldocchi 2014), the Bowen Ratio (Scott et al. 2006), sap flow, large aperture scintillometers, lysimeters and others (Allen et al. 2011). The eddy covariance method is considered the most reliable and state-of-the-art technique for accurate ET measurement at the field scale (Zhang et al. 2014; Saitta et al. 2020).

All methods mentioned above can cover small to medium-sized areas. The estimation of ET for larger areas is performed with the use of ET models. Several ET modeling approaches exist with different degrees of complexity (Chen & Liu 2020). The basis of those models is solving the energy balance equation using the well-known Penman-Monteith equation (Allen et al. 1998). ET models can be summarized into two different categories: (a) those that are based on the land surface temperature (LST), which is closely linked to ecosystem water status and thus ET (Bastiaanssen et al. 1998; Anderson et al. 2007; Mu et al. 2013b; Chen & Liu 2020) and (b) those that are based on the crop coefficient approach, the so-called FAO approach (Allen et al. 1998). The actual evapotranspiration was estimated as the product of reference evapotranspiration (ET0) multiplied by a factor called crop coefficient (Kc).

Both approaches have significantly benefited from the advantages that remote sensing is offering for large scale monitoring. The first category of models uses the thermal bands of Modis or Landsat satellites to estimate LST by applying quite sophisticated approaches (Anderson et al. 2012; Mu et al. 2013b). However, as satellite thermal bands are provided in lower spatial resolutions when the modeling applications refer to smaller areas, there is a need for downscaling with special algorithms (Agam et al. 2007; Gao et al. 2012).

The second category has been so far following a more semi-empirical approach. While reference evapotranspiration is estimated as a function of environmental parameters, the crop coefficients are empirical species, specific factors, fluctuating during the year according to the species' growth stage (Allen et al. 1998; Pereira et al. 2021). In the past few years, satellite vegetation indices have been recognized to offer great advantages in crop coefficient modeling (Pereira et al. 2021). So far, several studies have used satellite spectral indices to represent the seasonal fluctuation of crop coefficients, such as NDVI (Hunsaker et al. 2005; Drerup et al. 2017; French et al. 2020), or relative indices, like EVI and Soil Adjusted Vegetation Index (SAVI) (Guerschman et al. 2009; Glenn et al. 2011; He et al. 2019). Most of the studies that follow similar approaches refer to agricultural sites, which means that water availability is not a marginal factor, so NDVI and EVI can represent the phenological status. However, it is difficult to describe the water stress effects. In satellite primary productivity models that use the light use efficiency approach (Monteith 1972), the Land Surface Water Index (LSWI), (Gao 1996)) was widely used as a representative of ecosystem water status (Xiao 2004; Mahadevan et al. 2008; Chen et al. 2011; Bao et al. 2016). The index is related to the total canopy water content (Gao 1996), and thus it could potentially represent both the canopy growth stage and water stress. This approach can provide further benefits; the vegetation indices mentioned can be acquired with a high spatial resolution, such as Sentinel 2, making the applications more suitable for small or medium-sized areas.

Modeling applications have been mainly focused on agricultural sites or natural ecosystems. On the other hand, the applications in forest plantations on disturbed or degraded land, like black locust plantations (Kanzler et al. 2020), are generally missing.

Black locust (Robinia pseudoacacia L.) is a medium-sized deciduous species. Although it is considered an invasive species in Europe, it has already become a very common component of both economy and culture in many European countries (Vítková et al. 2017). Even though the use of the species in forest plantations has met a lot of criticism, the species has a long history of use in restoration projects due to its unique characteristics to adapt to marginal sites. It is a fast-growing, drought-tolerant and low demand species, able to grow in poor-nutrient areas (Nicolescu et al. 2020). Because of these unique characteristics, the species has become very common in reclamation projects for disturbed lands, such as former coal mines (Skousen & Zipper 2014; Carl et al. 2017; Vítková et al. 2017; Nicolescu et al. 2020). In Europe, black locust has been studied for a long time, with the studies being focused more on its economic value than ecology (Vítková et al. 2017). Over the last few years, many studies of the species have been mainly from the Loess Plateau in China (Wang et al. 2015; Jiao et al. 2018; Miyazawa et al. 2018; Zhao et al. 2018). Despite that, there are things about its ecology and ecophysiology that remain elusive.

The aim of this study is (a) to provide a medium-term presentation of measurements of daily ET fluxes for a tree plantation with the use of the eddy covariance method, (b) to assess the contribution of the main environmental driving parameters of the ET fluxes for the specific site and (c) to provide a simple, but accurate methodology for the extrapolation of the eddy tower measurements to greater scales with the use of a hybrid satellite and meteorological driven model, based on high spatial analysis vegetation indices from the Sentinel 2 satellite.

The site is located on a flat field of about 5 × 10⁵ m² (40.59 °N, 21.65 °E) that is part of a greater planted area of about 24 × 10⁶ m² in the restored areas of the Western Macedonia Lignite Center, in Northern Greece (Figure 1). The dominant species of the plantations is Robinia pseudoacacia, a deciduous broadleaf species, which accounts for 95% of the total area of the plantations. The plantations consist of mature and younger trees, due to the continuous restoration and plantation process of more than 20 years. The age of the trees at the study site is about 20 years. The mean tree height in the study site is 13.5 m and is relatively homogenous. Robinia pseudoacacia is a fast-growing species that can reach its maximum height before the age of 20 years (Nicolescu et al. 2020). The seasonal fluctuation of the Enhanced Vegetation Index (EVI) in the plantation, as estimated from the Sentinel 2 imagery, for a five-year period (2016–2020) is steady (see supplementary material). EVI is an index closely related to Leaf Area Index (LAI) (area of transpiration). Understory vegetation consists mainly of perennial grasses, with Cynodon dactylon being the dominant one. The growing period of the understory grass starts after leaf fall, as the leaves of black locust are very thin and delicate, resulting in a very quick decomposition, and it is active until the trees' leaf expansion at the beginning of May. The substrate consists of deposition material from the mining process that, in general, is very depleted and poor. It is characterized by alkaline pH and classified from moderately coarse to fine soils (Papadopoulos et al. 2015). Mean annual precipitation for the ten-year period in question (2010–2020) is 510 ± 156 mm, and mean annual temperature for the same period is 13.36 ± 0.92 °C, as estimated from the nearby meteorological station of the National Observatory of Athens (http://penteli.meteo.gr/stations/amyntaio/).

An eddy tower was installed at the end of July 2019, measuring continuously the exchanges of vapor, energy and CO₂ in the study area. The tower was installed in the middle of the study site in order to achieve a unidirectional footprint (no prevailing wind direction). The tower is equipped with an integrated IRGA analyzer and sonic anemometer (Irgason, Campbell Scientific Inc.) with an enhanced barometer (PTB 110, Vaisala Inc.), which perform CO₂/H₂O, three-component wind analysis and barometric pressure respectively and a full meteorological station, including, among others, a net radiometer (NR Lite2, Kipp and Zonen BV), a temperature/humidity probe (HygroClip2 Advanced, Rotronic Inc.) and three soil heat flux plates (HFP01, Hukseflux Thermal Sensors BV). The Irgason analyzer was placed in the top of the tower, at a total height of 17.5 m, providing us an average fetch radius of 200 meters. Soil water content (SWC) is measured at 15 cm depth. Raw data are collected from a datalogger (CR1000X, Campbell Scientific Inc.) at a frequency of 10 Hz and processed on-site via an integrated analysis software (EasyFlux DL, Campbell Inc.) and post-processed with EddyPro 7 (Licor GmBh) with the use of the standardized eddy covariance methodology (Aubinet et al. 1999). The default time step of data acquisition was 30 min. Daily values were obtained by integrating or averaging the instantaneous values, depending on the parameter. Data gap-filling and flux partitioning were done using a standardized methodology (Reichstein et al. 2005). Data analysis covered the period from August 2019 to September 2020.

All the above parameters are estimated by the eddy tower on half-hourly time steps and then scaled to daily time steps. By definition, air temperature and wind speed must be measured at the height of 2 meters or, if measured elsewhere, they have to be converted to that height. Our measurements were performed at the tower height (17.5 m) and no conversions were applied, to have a common reference with the other tower measurements.

For the calculation of vegetation indices, we used the free-distributed products of Copernicus Sentinel 2 (https://sentinel.esa.int/web/sentinel/missions/sentinel-2), which offer the advantage of high resolution (10 × 10 m) and frequent acquisition. For our analyses, we used L2A products. The average acquisition of cloud-free images during the study period for our study site was four days.

Using five bands in total (blue, green, red, nir and swir), we estimated the time series of the vegetation indices Land Surface Water Index (LSWI), Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI) and Normalized Difference Water Index (NDWI) (Table 1) for the whole study period.

All vegetation indices were estimated as the average ± standard deviation for the whole study site for cloud-free days, using the online analysis tool Sentinel-hub (https://www.sentinel-hub.com/). Sentinel 2 provides two SWIR bands, one at 1,565–1,655 nm and the other at 2,100–2,280 nm (https://sentinels.copernicus.eu). In this study, we used the SWIR1 (1,565–1,655 nm) band for the estimation of LSWI.

Due to the specific conditions in the study site; that is, the two distinctive states of the vegetation (seasonally dominated either by grass or black locust canopy), we used two separate crop coefficients for different periods of the year. For the period between the leaf fall and new leaf expansion (November to April), during which the site's vegetation was dominated by the understory grass, the value of Kc_max equal to 1.2 was used (Allen et al. 2011). For the period after leaf expansion until leaf fall (May-October), during which the vegetation was dominated by the canopy, the value of Kc_max equal to 0.75 was used, which was the average maximum value of the ecosystem's Kc during the period of full leaf expansion, estimated as the ratio of ET/ET0 during the corresponding period, according to FAO methodology (Allen et al. 1998; Pereira et al. 2021). The value 0.75 for Kc_max used is inside the range of values for Kc of deciduous trees (0.45–1.2) proposed by FAO (Allen et al. 1998; Pereira et al. 2021). For the intermediate period from the beginning of leaf expansion until full development (April-May), during which both understory vegetation and canopy are active, the average of the two values was used. The vegetation index used as a scalar for Kc_max (Equation (3)) was better correlated to the measured ET fluxes (see results). As min and max values in Equation (3) minimum and maximum values of the indices during the study period were used and as the running values, the calculated indices, interpolated for the interval periods.

The second approach was similar to the one described above, with the difference that instead of ET0 we used the H₂O equivalent for net radiation (mm d⁻¹) as input to Equation (2). ET0 is a quantity that depends on many environmental parameters, as extensively described before. All the parameters do not have the same impact on actual ET. The subtraction of soil heat flux from net radiation should ideally be equal to the sum of latent and sensible heat flux (LE + H) which determines the total evaporative demand. We tried to further minimize the input parameters. An alternative version of the model was utilized. Total evaporative demand is not represented by reference evapotranspiration, but by the water equivalent of net radiation (mm d⁻¹). However, as the sum of LE + H are usually less than net radiation – global radiation (Rn-G) (Allen et al. 2011) and G is small compared to Rn, the water equivalent of Rn was multiplied by 0.88, a factor that resulted from the solving of the energy balance equation for our study site (see supplementary material).

The third approach was the derivation of an empirical formula from the regression of the measured daily evapotranspiration values with the satellite vegetation indices mentioned in the previous section (Table 1). Although this is the simplest approach, it lacks a mechanistic interpretation. While it can be useful for the specific site, it would be quite uncertain to use it in transferable applications.

The assessment of the importance of relative contribution of each independent variable (net radiation (Rn), wind speed (WS), temperature (T) and vapor pressure deficit (VPD) on ecosystem evapotranspiration was done by means of P values and standardized beta coefficients after the application of multiple regression analysis.

The seasonal fluctuation of meteorological data from the eddy tower is presented in Figure 2 and was used as inputs in Equation (1) to estimate reference evapotranspiration. The peak of net radiation (Rn) was observed during the period between mid-June to mid-July. Average Rn values were slightly higher during August 2020 compared to August 2019. Differences between 8/2019 and 8/2020 were observed at air temperature (T) and vapor pressure deficit (VPD) (Figure 2(b) and 2(c)). Both quantities were lower during 8/2020 compared to 8/2019 by about 2 °C and 0.3 kPa respectively. The precipitation and SWC are presented in Figure 2(d). The total precipitation amount for the whole study period (1/8/2019–27/9/2020) was 622 mm. Summer precipitation for 2020 (June – August) was quite high (170.8 mm), which accounted for 27.5% of total precipitation for the whole study period. SWC was kept high from the end of autumn until the middle of spring, when a rapid decrease was observed from the end of April until the end of May. In general, summer SWC values were quite high, as well as the precipitation. However, interannual differences were also observed, as the average % SWC during 8/2019 was 9.09 ± 1.55, while in 8/2020 23.54 ± 4.99 and the corresponding precipitation values were 3.1 and 36.9 mm, respectively. Finally, wind speed values did not appear to follow a clear seasonal pattern or have great differentiations, except for some higher values during the winter (data not shown).

The seasonal fluctuation of both reference (ET0) as estimated by Equation (1) and actual evapotranspiration (ET), measured by the tower, is presented in Figure 3. The overall seasonal fluctuation of both quantities is generally following the annual pattern of net radiation and SWC (Figure 2). ET0 value is higher than the ET during the whole measured period. The lower annual rates are met during January, while two peaks can be found, one from May to the beginning of June and the other from July until the middle of August. Between those periods there is an intermediate period of a sharp decline in the activity, which reflects the significant drop in rainfall and the decrease in SWC. The ET rates are considerably high from February to April, before the leaf expansion, which can be attributed to the well-developed understory grass.

The period of the first peak (May-June) is similar to the period of ET maximization for the same species in the Loess Plateau in China (Miyazawa et al. 2018). The trees of the referred study were grown on an arid to semi-arid region with a different precipitation pattern and that's probably the reason for this very narrow peak. The peak during the beginning of May in our study site could be considered quite early, as during May the canopy was still under development and full leaf expansion is met about the middle of June (data not shown). It seems this species and the whole ecosystem, consequently, follow a strategy of very luxury water consumption during the canopy development. The result of this is rapid soil drying until the beginning of June, despite the considerable amount of precipitation during that period (Figure 2(d)). The excessive water consumption of black locust and the resulting decline in soil moisture has been reported in another study from a relative ecosystem in China (Liang et al. 2018). The estimated weighted annual sum of ET with eddy covariance for Robinia pseudoacacia is 595 mm, a value that is very close to the respective total sum of precipitation. Although this value is comparable to many deciduous broadleaf sites of the fluxnet network (fluxnet.org) (Fatichi & Pappas 2017), there are not many sites with the same species that are comparable, as was mentioned previously. In a relevant study, the average total ET of a two-year study for R. pseudoacacia was estimated to be about 208 mm y⁻¹. However, the average amount of precipitation was quite lower, about 262 mm y⁻¹ (Jiao et al. 2018). Another study, (Fu et al. 2021) estimated the annual evapotranspiration for an R. pseudoacacia plantation under non-water stress conditions at 576.88 and 848.55 mm y⁻¹ for a dry and a wet year respectively.

The assumptions mentioned before are made with a quite high degree of uncertainty, since our results refer to the total ecosystem fluxes and any partitioning of ecosystem ET to transpiration and soil evaporation was not applied. Although generally the contribution of transpiration to total ecosystem evapotranspiration in closed canopies is estimated to be up to 90% (Jasechko et al. 2013; Good et al. 2014; Li et al. 2019), in a relevant study in a R. pseudoacacia plantation, the corresponding proportion was estimated between 9.7 and 28.8% (Jiao et al. 2018), which can be considered quite low.

Measured ET values are strongly correlated with all studied vegetation indices, with coefficients of determination ranging from 0.51 to 0.77. The strongest correlated index was the LSWI (R² = 0.77, Figure 4(a)). LSWI was used as the input parameter to describe the Kc scalar (Equation (3)) for the models based on reference evapotranspiration (Equation (2)). However, as the initial idea was to use a scalar for Kc that would be representative of both phenological and water status, the performance of LSWI versus ET was investigated under different water status scenarios. Given that daily ET seems to co-fluctuate seasonally with SWC (Figures 2 and 3), the effect of SWC to ET was assessed by comparing the respective values for the period of full leaf expansion (June – August, Figure 5(a)). This period limitation was followed to avoid the misinterpretation or noise that would come from data from the non-growing period. According to our findings, SWC is related to daily ET via an exponential function (Figure 5(a)) with a very high coefficient of determination (R² = 0.76, P < 0.001). Based on this relationship, an apparent threshold of 15% SWC can be identified, below which ET values are strongly increasing with the increase of SWC and reach a plateau at about the value of 15% SWC. Following that threshold, the LSWI – ET values were separated into two sub-datasets, one that corresponded to SWC values greater than 15% and the other to SWC values less than 15% for the whole period of study, including the non-growing period. Regression analysis with the use of an exponential function was performed for the two datasets independently, as well as for the whole dataset together (Figure 5(b)). Based on the results of this analysis, ET is related to LSWI with exponential functions. Although ET values corresponding to SWC > 15% are significantly higher than the relative values that correspond to SWC < 15% for the same LSWI range, the use of a single exponential function applied to the whole dataset can be considered adequate to describe the ET dependence on LSWI for the whole SWC range, as nearly all values are inside the corresponding 95% prediction band (Figure 5(b)). Consequently, we can assume that LSWI is an appropriate index for the description of both the phenological stage and water status of our studied ecosystem. The accurate representation of water availability effects on ecosystems performance is one of the most important, yet demanding and uncertain procedures. Our findings highlight the great potential of LSWI to be used as a proxy for water limitations in our studied ecosystem and could act as a driver for further implications in other ecosystems.

Pearson R and regression coefficients R² between ET and ET0 daily values and the corresponding environmental parameters and LSWI are presented in Table 2. ET and ET0 are very strongly correlated with Rn, VPD, T and LSWI, while the correlation with wind speed is quite weaker for both quantities. The latter can be explained from the relative constant wind speed conditions throughout the study period. ET is better correlated with net radiation and LSWI, rather than ET0.

Both ET and ET0 are simultaneously affected by all the independent parameters mentioned above. The one-to-one comparison can lead to false conclusions about the real effect of each parameter on the final estimates. For example, net radiation (Rn) and temperature (T) follow similar seasonal fluctuation patterns, so their true relative impact cannot be distinguished in that way. Additional multiple regression analysis was applied, in order to estimate both the significance and the relative effect of each studied environmental parameter to ET and ET0. This was assessed by the P values and the standardized beta coefficients for each parameter. The results of the analysis are presented in Table 3.

According to the multiple regression analysis results, ET0 was significantly affected by all environmental parameters (P < 0.001), while the parameters that have greater effect on ET0 are net radiation and VPD, with beta values 0.64 and 0.39, respectively. The effect of temperature and wind speed are quite weaker. All independent parameters were positively correlated to ET0.

For actual ET, when LSWI was not included in the analysis, all environmental parameters were significant (P < 0.001) and Rn was the parameter with the strongest effect on ET. The second most important parameter was T. VPD was also found to have a significant effect; however, contrary to ET0 on this occasion, the effect is negative. The results of the analysis were altered when LSWI was included as an independent parameter in the analysis. In this case, both air temperature and VPD had no significant correlation to ET, as their P values are 0.27 and 0.71, respectively. The only important parameters are Rn, WS and LSWI. The parameter with the strongest effect was LSWI, followed by net radiation and wind speed.

Regarding Rn, these results are somehow expected, as Rn was the strongest component in the energy balance equation and thus determined the total evaporative demand. However, when T and VPD are analyzed, things are more complicated. Firstly, VPD is dependent on air temperature, so the use of both parameters in the same analysis can create some uncertainty. Secondly, while the increase of VPD increases the total evaporative demand, its effect on actual evapotranspiration could be either positive or negative (Massmann et al. 2019) as the increase of VPD leads to decrease in stomatal conductance (Dewar 2002).

According to the above, the model can be further simplified by using the water equivalent of net radiation (Rn) instead of reference evapotranspiration as a proxy for evaporative demand. This led to a model that relies on simplicity, without compromises in accuracy (R² = 0.89, RMSE = 0.53, NSE = 0.80). This can provide further advantages, as net radiation is a parameter that can be easily measured in the field or estimated from relevant satellite maps. The problem in this approach is the energy closure. The sum of LE + H does not equal to Rn-G, with an uncertainty that can be as high as 30%, with usual values around 10–20% (Wilson et al. 2002; Allen et al. 2011; Burba 2013). In the case of our studied ecosystem, the relative value is estimated around 12% (see supplementary material).

Minimizing the input parameters, without making great compromises in accuracy, is of great importance in ecosystem modeling. According to our findings, not all environmental parameters have the same impact on ET fluxes for the studied ecosystem. Net radiation and LSWI can be adequate in order to provide an adequate estimation of ET on seasonal basis. Considering that net radiation is a quite easily measured or estimated parameter, that satellite imagery offers many advantages in its estimation (Hou et al. 2014; Verma et al. 2016; Chen et al. 2020). Its use along with LSWI could meet further applications in evapotranspiration modeling (Barr et al. 2012).

The seasonal fluctuation of measured against modeled ΕΤ values, according to the three different modeling approaches, as well as the comparison of the measured and modeled values, is presented in Figure 6. The modeled values derived from Equation (2), using ET0 as a proxy for evaporative demand, are presented as Model 1, while the corresponding values that use the water equivalent of net radiation as a proxy for evaporative demand are presented as Model 2. Finally, the modeled ET values that come from the empirical regression of measured ET with LSWI, are presented as Model 3. All the three models show very good seasonal co-fluctuation and agreement with the measured values for their 8-day running averages (Figure 6(a)–6(c)), with comparable R² values, ranging from 0.85 to 0.89 (Figure 6(d)–6(f)). The same applies when RMSE and NSE are considered, with values ranging from 0.45–0.53 and 0.80–0.85, respectively. It must be mentioned that typically this comparison can be considered valid only for the cases of Model 1 and Model 2, because the regression analysis in the case of Model 3 comes from the same dataset with which the model was evaluated, so typically that can be considered as autocorrelation. However, our scope in this case is to present how the modeled values of the empirical function co-fluctuate seasonally with the measured values and what are the advantages or disadvantages compared to the more mechanistic models.

In general, ET models have been reported to perform quite satisfactory in several studies with R² coefficients ranging from 0.3 to 0.86 for crops (Battude et al. 2017; He et al. 2019; French et al. 2020) and 0.12–0.90 in natural forests (Nagler et al. 2005; Verstraeten et al. 2005; Cleugh et al. 2007; Velpuri et al. 2013; Yao et al. 2017). According to that, our modeling approach resulted in a model performance among the higher range.

However, the advantages of the developed model are not mainly concerned with its performance, but its simplicity. Although several models have used similar approaches; that is, to use satellite indices to model the crop coefficient, the most commonly used indices were NDVI and EVI (Hunsaker et al. 2005; Guerschman et al. 2009; Glenn et al. 2011; Drerup et al. 2017; He et al. 2019; French et al. 2020), which can adequately describe the phenological stage of the vegetation but are poorer for its water status. On the other hand, LSWI is related to the total vegetation water content, and as a result to both phenological stage and water status (Gao 1996). Additionally, LSWI has been found to have a strong correlation with LST (Mokhtari et al. 2019; Alexander 2020) and thus can be used as a substitute of thermal spectrum (Sadeghi et al. 2015). Therefore, the use of LSWI, as a scalar to crop coefficient, could potentially work as a bridge between LST and crop coefficient based models. This is even more important if we considered that the studied ecosystem is not an agricultural site or a natural forest, but a plantation grown on a degraded site.

In this study we presented the results of more than one-year measurements and modeling of ecosystem evapotranspiration for a Robinia pseudoacacia plantation grown after restoration of highly disturbed land.

Reference evapotranspiration estimated from Equation (1) (ET0) has given a higher value from actually measured evapotranspiration (ET) from an eddy tower during the whole study period. A high value of ET was recorded in May to middle August. SWC is related to daily ET via an exponential function. The relationship between satellite spectral indices and the corresponding daily evapotranspiration values showed that LSWI was the most higher correlated index. From the environmental parameters, LSWI and Rn were those which had strong effect on ET. ET is dependent on LSWI and it is an appropriate index for the description of both phenological stage and water status of the studied ecosystem.

The results of this study contribute to the better understanding of the behavior of Robinia pseudoacacia, as well as the ecosystem processes of forest plantations grown on disturbed sites, and potentially they can lead to further investigation of the potential of minimization of the input parameters on ET models, without great compromises on model performance and accuracy.

In general, although R. pseudοacacia is already a species with an extensive distribution in Europe and worldwide, there are still many things to be studied concerning its ecology and behavior, and its contribution to the global carbon and water cycle. We believe that our findings can have a significant contribution in this direction.

The authors kindly acknowledge the Hellenic PPC (Public Power Corporation) S.A. for the provision of assisting personnel and necessary infrastructure during field campaigns and data collection. Special thanks are due to Lamprini Patmanidou and Simela Andreadi and their teams.

This research has been co-financed by the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship, and Innovation, under the call RESEARCH–CREATE–INNOVATE (project: COFORMIT. T1EDK-02521).

Data cannot be made publicly available; readers should contact the corresponding author for details.