Improving remote estimation of winter crops gross ecosystem production by inclusion of leaf area index in a spectral model

The hysteresis of the seasonal relationships between vegetation indices (VIs) and gross ecosystem production (GEP) results in differences between these relationships during vegetative and reproductive phases of plant development cycle and may limit their applicability for estimation of croplands productivity over the entire season. To mitigate this problem and to increase the accuracy of remote sensing-based models for GEP estimation we developed a simple empirical model where greenness-related VIs are multiplied by the leaf area index (LAI). The product of this multiplication has the same seasonality as GEP, and specifically for vegetative periods of winter crops, it allowed the accuracy of GEP estimations to increase and resulted in a significant reduction of the hysteresis of VIs vs. GEP. Our objective was to test the multiyear relationships between VIs and daily GEP in order to develop more general models maintaining reliable performance when applied to years characterized by different climatic conditions. The general model parametrized with NDVI and LAI product allowed to estimate daily GEP of winter and spring crops with an error smaller than 14%, and the rate of GEP over- (for spring barley) or underestimation (for winter crops and potato) was smaller than 25%. The proposed approach may increase the accuracy of crop productivity estimation when greenness VIs are saturating early in the growing season.


INTRODUCTION
Leaf area index (LAI) as well as parameters describing carbon dioxide (CO 2 ) exchange between plants and the atmosphere such as net ecosystem production (NEP) and gross ecosystem production (GEP) are key biophysical parameters, which are commonly applied to qualitatively and quantitatively characterize the status of vegetation canopies. Numerous studies have confirmed that based on LAI the intensity of photosynthesis, transpiration, and productivity of plants might be assessed (Borge & Leblanc, 2001;Breda, 2003;Zarco-Tejada, Ustin & Whiting, 2005). Hence, LAI can be used as a proxy of plant growth, biomass and yield, as well as carbon dioxide fluxes exchanged between the ecosystem and the atmosphere (Breda, 2003;Glenn et al., 2008). GEP reflects the total amount of CO 2 assimilated by plants in photosynthesis (Waring, Landsberg & Williams, 1998). It depends on the amount, type and physiological condition of plants, but also on climate and habitat conditions (Baldocchi et al., 2001;Keyser et al., 2000). GEP flux analysis is of importance of studies concerning carbon assimilation efficiency at leaf, plant and ecosystem levels. However, due to the limitations of measurement methods, GEP cannot be directly measured in situ. State-of-the-art eddy covariance (EC) systems installed on flux towers make it possible to measure only the net fluxes of CO 2 (NEP) (Urbaniak et al., 2016). NEP is defined as a balance of the processes of CO 2 exchange between the ecosystem and the atmosphere and is expressed as a difference between GEP and the total amount of CO 2 released by the ecosystem to the atmosphere (ecosystem respiration, R eco ) (Law et al., 2002).
Remote sensing studies of global vegetation phenology started in 1979 when the meteorological satellite data of Advanced Very High Resolution Radiometer (AVHRR) became available (Goward & Huemmrich, 1992). Since then, the following satellite missions (MODIS, LANDSAT, Sentinel-2) have been providing data which allow remote sensingbased estimation of LAI, fraction of PAR absorbed by plants (fAPAR) and GEP of biomes and ecosystems across the globe with higher spatial and temporal resolutions as well as higher accuracy (e.g., Running et al., 2004;Gao et al., 2017). Traditionally, GEP is estimated as a function of vegetation indices (VI ) related to canopy greenness or based on models including in their formulation also Light Use Efficiency (LUE) and/or Photosynthetic Active Radiation (PAR) terms (e.g., Gitelson et al., 2012;Rossini et al., 2012;Sakowska et al., 2014). Many studies have highlighted that simple greenness-related VIs can be successfully used for remote sensing-based estimation of LAI (e.g., Gitelson et al., 2003), chlorophyll content (e.g., Gitelson et al., 2005), fAPAR (Sims et al., 2006), fractional vegetation cover (e.g., Glenn et al., 2008) and GEP (e.g., Prince & Goward, 1995). Among these greenness indices, Normalized Difference Vegetation Index (NDVI ) has been the most commonly applied, even though it tends to saturate under conditions of moderate-to high aboveground biomass (e.g., Gitelson, 2004). For this reason, a big effort was undertaken to develop new NDVI -type indices that would not only minimize the soil background influences (e.g., Soil Adjusted Vegetation Index, SAVI, Huete, 1988), but would also overcome the saturation problem in biophysical parameters estimation (Gitelson, 2004). According to existing studies, Modified Simple Ratio (MSR, Chen & Cihlar, 1996), Renormalized Difference Vegetation Index (RDVI, Roujean & Breon, 1995), Wide Dynamic Range Vegetation Index (WDRVI, Gitelson, 2004) or Enhanced Vegetation Index (EVI Huete et al., 2002;Rahman et al., 2005) are more linearly related to fAPAR, LAI, or GEP and they allow to estimate these biophysical parameters with higher accuracy. Due to the limited sensitivity of ''greenness'' indices to a short-term stress which may not impact the chlorophyll content, the Photochemical Reflectance Index (PRI) was introduced (Gamon, Penuelas & Field, 1992). PRI may be an indicator of the LUE in the process of photosynthesis (Gamon, Penuelas & Field, 1992;Goerner et al., 2011;Penuelas, Filella & Gamon, 1995) and has been used in GEP and LUE estimations at leaf, plant and ecosystem levels (e.g., Rossini et al., 2012;Cheng et al., 2014;Gitelson, Gamon & Solovchenko, 2017).
In order to obtain remote sensing-based model which allows to estimate daily GEP of crops independently from the type of the crop and climatic conditions with reliable performance, we studied a 3-year dataset (2011-2013) consisting of spectral and biophysical data for two winter (wheat and rye) and two spring (barley, potato) crops. Our specific objectives were to test (1) the accuracy of daily GEP estimations with remote sensing-based models fed with different VIs and VIs*PAR and a model based on LAI, and (2) whether the accuracy of GEP estimations increases when products of VI and LAI are included in the model. Besides, considering that for some crops the hysteresis of the relationships between VIs and biophysical parameters is observed (Peng et al., 2017), we aimed at developing a simple linear empirical model based on VIs and LAI in order to reduce hysteresis of the VIs vs. GEP relationships between vegetative and reproductive phases of crop development cycle and to increase the accuracy of daily GEP (GEP d ) estimations of croplands. According to our knowledge this is also the first study in which CO 2 fluxes measured with chambers are combined with both LAI and VIs.

Experimental site
Measurements were conducted at the Brody Experimental Station (52 • 26 N, 16 • 18 E) on plots of the long-term experiment that has been conducted since 1957 by the Department of Agronomy, Poznań University of Life Sciences, Poland (Blecharczyk et al., 2016). Crops were grown in the crop rotation and monoculture systems under 11 different fertilization regimes (no fertilization, manure, manure + NPK, NPK + Ca-CaO, NPK, NP, NK, PK, N, P, K). The measurements presented in this paper were performed on four crop species: potato (var. Wineta), spring barley (var. Nadek), winter wheat (var. Turkis) and winter rye (var. Dankowskie Zlote), grown in a seven-year rotation (potato → spring barley → winter triticale →1-and 2-year alfalfa → winter wheat → winter rye). The investigated crops were fertilized with NPK (90 kg N ha −1 a −1 , 60 kg P 2 O 5 ha −1 a −1 , 120 kg K 2 O ha −1 a −1 ) with an addition of Ca-CaO (1.5 Mg CaO ha −1 a −1 ) and grown in 6 × 11 m plots separated by 0.5 m wide bare soil stripe. The annual mean air temperature of the study area is 7.9 • C, while the annual precipitation sum is 571 mm (average for 1959-1999). The soils are classified as Albic Luvisols developed on loamy sands overlying loamy material (Majchrzak et al., 2016).

Chamber CO 2 fluxes and flux modelling
Measurements of CO 2 fluxes (NEP, R eco ) were taken using a closed dynamic portable chamber system, which consisted of a transparent and non-transparent chambers to measure NEP and R eco fluxes, respectively (Chojnicki et al., 2010;Juszczak et al., 2013). Measurements were made in two subplots for each crop. The transparent chamber was made of three mm-thick Plexiglas (Evonik Industries, Darmstadt, Germany), as this material has a high solar radiation transmittance (approximately 90%, Acosta et al., 2017;Hoffman et al., 2015). The non-transparent chamber was made of three mm-thick white PVC to ensure dark conditions inside the chamber. The chambers had dimensions of 0.78 × 0.78 × 0.50 m and a total volume of 0.296 m 3 . In case of winter rye and winter wheat, the extensions of 0.5 m height were used (made of the same material) in order to adopt the height of chamber to the height of the canopy. During the measurements, chambers were placed on square PVC collars (0.75 × 0.75 m), inserted into the soil just after sowing. The insertion depth of the collars was 15 cm. The chambers were equipped with a set of computer fans (1.4 W; 1,500 rpm each) mixing the air and a vent to equilibrate pressure in the chamber headspace. The air temperature inside the chamber headspace was measured with a radiation-shielded thermistor (T-107; Campbell Scientific, Logan, UT, USA) at a height of 0.3 m, or 0.8 m for short and tall chambers, respectively (Juszczak, Acosta & Olejnik, 2012;Juszczak et al., 2013). The CO 2 concentration changes in the chamber was measured using LI-820 gas analyzer (LI-COR Inc., Lincoln, NE, USA). The air was circulated between the chamber and the analyzer in a closed loop with the flow rate of 0.7 l min −1 . In order to keep the air temperature inside the chamber headspace stable, the transparent chamber was cooled with a passive system described in Acosta et al. (2017).
Chamber measurements were taken every 3-5 weeks throughout the entire year (including winter) on cloudless days from sunrise to late afternoon (Uździcka et al., 2017). However, when any clouds appeared (usually in the afternoon) a particular attention was paid to perform measurements at stable PAR conditions. Overall, 37 chamber campaigns were conducted in the years 2011-2013. Measurements of NEP and R eco were taken on each of the soil frames several times per day (from five to 12, depending on the daytime length). A single NEP measurement took 120 seconds, and the subsequent NEP measurement at the same plot was taken before the value of incoming PAR changed by more than 150 µmol m −2 s −1 . A single measurement of R eco took 180 seconds and the succeeding measurements were taken before the soil temperature at the same plot changed by more than 0.5 • C. The CO 2 flux (F ) in µmols m −2 per time unit (t ) was calculated from the gas concentration change in the chamber headspace ( C t ), the chamber volume (V ) and the enclosed soil area (A) from the following equation (Urbaniak et al., 2016): where M v (m 3 mol −1 ) is the molar volume of air at a given chamber air temperature and pressure. The determination coefficient (r 2 ) was calculated for each single chamber closure time and if r 2 < 0.8, the fluxes were excluded from the analyzes. To ensure that good quality near-zero fluxes were not erroneously excluded by this criteria, all rejected fluxes were visually inspected. Fluxes of NEP, R eco and GEP for measurement days and periods between campaigns were calculated using a simple empirical model described by Drösler (2005) and farther elaborated by Hoffman et al. (2015). For this purpose, for each day of measurements, the relationships between measured R eco and temperature were established by fitting to the campaign-specific flux dataset the temperature dependent Arrhenius-type respiration model of Lloyd & Taylor (1994). Using the parameters of this model, R eco values at the time of the NEP measurements were estimated based on the measured temperatures (Juszczak et al., 2013). In the next step, based on such estimated R eco and measured NEP, GEP was calculated according to the formula GEP = NEP + R eco . Subsequently, the calculated GEP fluxes were correlated with measured PAR, fitting to the campaign-specific GEP dataset a rectangular hyperbolic light response Michaelis-Menten kinetic model (Michaelis & Menten, 1913). The R eco and GEP model parameters were interpolated linearly for the periods between the campaigns with a 30-minute step, so that, based on the continuous time series of measured PAR and temperature (means for 30-minute periods), GEP and R eco was calculated. NEP was calculated from the formula NEP = GEP − R eco . Daily sums of GEP (so called daily GEP (GEP d )) were calculated as a sum of all 30-minute GEP fluxes estimated for each day with the flux model in between sunrise to sunset.

Measurements of LAI and spectral characteristics of the crops
LAI and multispectral data were collected during the growing season, from March to October, at one-to two-week intervals. At the spring barley plots measurements started in the middle of April (after sowing), whereas in case of potatoes in the 2nd week of May Gamon, Penuelas & Field (1992) WDRVI WDRVI = α * ρ850−ρ670 α * ρ850−ρ670 Gitelson (2004) Notes. ρ, reflectance at a given wavelength; NDVI, Normalized Difference Vegetation Index; SAVI, Soil Adjusted Vegetation Index; PRI, Photochemical Reflectance Index; WDRVI, Wide Dynamic Range Vegetation Index.
(after planting). Only the data collected in the period between April and the 2nd week of August (just after harvesting) was considered in the analyzes, hence we did not present nor analyzed data collected after sowing of winter crops and after harvest. The dates of the measurement campaigns were selected so that these measurements could overlap with chamber measurements of CO 2 exchange. However, when the weather conditions were not stable due to appearing clouds, the measurements were repeated on the first sunny day following the chamber measurements. 36 measurement campaigns were organized in the 3-year period (2011)(2012)(2013). Spectral measurements were carried out only on sunny days, always around the solar noon (between 10:00-14:00). Measurements of LAI and reflectance were taken at each plot in three replications, always in the same locations. LAI was measured by means of the SunScan system (Delta-T Devices, Cambridge, UK). The spectral characteristics of the surface of plant-covered plots were measured using two 4-channel SKR1850 sensors (SKYE Instruments Ltd., Llandrindod Wells, UK) mounted on a portable SKL908 device (Spectrosense2+). Incident and reflected radiation was recorded at central wavelengths of 531, 570, 670 and 850 nm with 10 nm bandwidths. Next, applying the methodology developed by SKYE instruments Ltd. (SpectroSense2+ Manual), vegetation indices (NDVI, SAVI and PRI ) were calculated (Table 1) according to the formula: where VI is a vegetation index, Z is a ratio sensitivity of reflected R1 r :R2 r ; X and Y are incident readings for R1 i and R2 i , respectively (in µmol m −2 s −1 ), while R1 r and R2 r correspond to reflected signal readings for wavelengths R1 and R2 (in nanoamps). For NDVI and SAVI, R1 and R2 correspond to 850 nm and 670 nm wavelengths respectively, while for PRI they correspond to wavelengths 570 nm and 531 nm. To calculate SAVI, the above formula was modified according to equation provided in Table 1. Wide Dynamic Range Vegetation Index was calculated from NDVI from the equation: Gitelson, 2005). In order to express values of this VI in positive numbers, sWDRVI was calculated from equation (sWDRVI = (WDRVI + 1)/2). Table 2 The models tested for GEP d estimation in the present study.

Models for GEP estimations
Daily GEP (GEP d ) were estimated based on linear regressions assuming a direct linear relationship between GEP d and LAI (model 1), GEP d and VIs (model 2), GEP d and a product of VIs and mean daily PAR calculated for the time between the sunrise and sunset −PAR d (model 3), as well as GEP d and a product of VIs and LAI (model 4) ( Table 2). All the models were tested based on the combined multi-year (2011-2013) dataset for each crop species separately, as well as for all cereals: winter wheat, winter rye and spring barley (i) and cereals and potatoes (ii) considered together in order to develop the general models for GEP d estimations for croplands. The general models with the best goodness of fit were then tested for each crop independently.

Statistical analysis
Pearson's correlation analysis was used to test the significance of the relationships between GEP d and (i) LAI ; (ii) VIs (NDVI, SAVI, WDRVI, PRI ); (iii) VIs*PAR d ; and (iv) VIs* LAI. Analyzes were conducted for each crop independently and for the combined datasets of cereals, as well as for cereals and potatoes considered together. Each of the four analyzed models' coefficients were found by fitting each model against GEP d . Goodness of fit statistics (coefficient of determination, R 2 ; root means square error, RMSE in gCO 2 -C m −2 d −1 ; and normalized root mean square error, NRMSE in %) were computed to compare the performance of the models.
In order to determine if there are significant differences in relationships between GEP d and VIs between the vegetative and reproductive phases of crop development, the two-sample t -test approach was applied. The differences between analyzed relationships were considered to be significant if p-value obtained from the test was lower than 0.05.
Due to limited number of data (n equals from 21 to 26, depending on the crop and VI ), validation of the best performing general models developed based on the combined dataset for cereals and potatoes was conducted for each crop species separately by comparing GEP d retrieved from chamber CO 2 fluxes with GEP d estimated based on the spectral model individually for each crop.

RESULTS
The seasonal variations of meteorological conditions for the analyzed growing seasons are presented in Fig (Table 3).
Maximum values of NDVI and SAVI were observed nearly at the same time as the peaks of LAI and GEP d (data shown in File S2). Moreover they were also lower in 2011 than in the other two years, even though the observed differences in maximum VIs values were less prominent than the differences in the analyzed biophysical parameters (Table 3). All these data clearly indicate the effect of drought which occurred in the late spring -early summer of 2011, when sums of precipitation were 50% smaller than the precipitation observed in the same period of 2012 and 2013.
The analysis of linear regression revealed that LAI explained minimum 60% of the variability in GEP d (Table 4). NDVI and SAVI explained from 52% to 72% of the variability in GEP d for winter crops and up to 81%-91% for spring crops, when crops were considered separately. For crop-combined dataset, which consisted of the data of all the analyzed crops, NDVI and SAVI explained 50% to 65% of the variability in GEP d . SAVI -based models worked better only in case of models developed for winter rye, while in case of other crops it did not lead to more accurate estimations of GEP d . Moreover, the SAVI -based models (R 2 = 0.50, NRMSE = 18.24%) developed for all the crops together were less accurate than NDVI-based models (R 2 = 0.65, NRMSE = 15.29%). The sWDRVI, which was expected to be more linearly correlated with GEP d , explained between 50% (winter rye) to 86% (spring barley) of the GEP d variability and did not improve the accuracy of GEP d estimations in the crop-combined model (R 2 = 0.65, NRMSE = 15.60%). The highest accuracy of model 2 was obtained for potato (if based on NDVI ; R 2 = 0.91, NRMSE = 10.94%). Inclusion of PAR d into model 3 did not improve estimations of GEP d for any of the investigated crops nor for the crop-combined datasets.
The inclusion of LAI into the ''VIs''-based models resulted in a general increase of their performance in case of the winter crops (Table 4). The highest increase of the accuracy of GEP d estimations was found for NDVI -based models. For winter wheat, RMSE decreased from 4.31 to 3.56 gCO 2 -C m −2 d −1 , while NRMSE decreased from 19.83% to 16.38% after inclusion of LAI into NDVI -based model. For winter rye, the changes were even more prominent -RMSE decreased from 4.72 to 3.14 gCO 2 -C m −2 d −1 , whereas NRMSE decreased from 20.15% to 13.42%. Inclusion of LAI into NDVI -based models developed for spring barley and potatoes led to a decrease of the accuracy of GEP d estimations (Table  4). For spring barley RMSE and NRMSE increased from 1.73 to 1.79 gCO 2 -C m −2 d −1 and from 12.33% to 12.70%, respectively. The highest reduction of the model accuracy was observed for potatoes -RMSE and NRMSE increased from 1.26 to 2.10 gCO 2 -C m −2 d −1 and from 10.94% to 18.30%, respectively.
For more general models developed for the cereals and crop-combined datasets the inclusion of LAI into the VIs-based models also resulted in an improvement of model performance. RMSE and NRMSE of NDVI and LAI -based models were the smallest and decreased from 4.09 to 3.42 gCO 2 -C m −2 d −1 and from 17.02% to 14.39% for cerealscombined models, and from 3.67 to 3.26 gCO 2 -C m −2 d −1 and from 15.29% to 13.57% for crop-combined datasets, respectively (Table 4). That is why we used the most general and the best fitting NDVI*LAI model developed for all the crops together (cereals & potatoes, Fig. 3) to estimate GEP d values for each crop separately (Fig. 4). Although there is a good agreement between observed and predicted GEP d values for all the crops (Fig. 3), GEP d estimated for winter crops and potatoes was underestimated, while GEP d of spring barley was overestimated, but the rate of over-or under-estimation did not exceed 25% (Fig. 4).

DISCUSSION
Most of the remote sensing-based models to estimate GEP of croplands are rather cropspecific. They were developed for maize, soybean (Gitelson et al., 2012), rice (Inoue et al., 2008), wheat (Wu et al., 2009), rye, barley, potato (Uździcka et al., 2017 and in the majority of studies it has not been tested if these models can be directly applied (without reparametrization) for estimation of CO 2 uptake for other crops. Here we presented a more general and robust approach based on combining the datasets for the two winter and two spring crops together and we proved that the accuracy of crop-combined models is not different from those developed for winter crops on a separate basis (NRMSE of the best crop-combined models are in range of 13-16%), but it may be lowered compared to simple VI -based models developed for spring barley and potato (NRMSE is in range of about 11% for spring crops, whereas NRMSE of the best crop-combined model is 13.57%, Table 4). The weakness of our study is related to the applied method for spectral properties measurements. The Spectrosense 2+ measuring system with 4-channel upward-and downward-looking multispectral radiometers allows to measure incident and reflected radiation only in four most commonly used bands (NIR, RED and GREEN wavelengths). Due to missing measurements at red-edge, blue and other specific wavelengths we were not able to calculate many other important indices more sensitive to medium to high biomass (e.g., NDVI red-edge (Gitelson & Merzlyak, 1994), Normalized Difference Structural Index, NDSI (Vescovo et al., 2012)) or EVI (Huete et al., 2002), which could help to overcome the saturation effect of classical greenness indices. Besides NDVI, which tends to saturate asymptotically under moderate-to-high biomass conditions (Huete et al., 2002;Gitelson et al., 2003), we analyzed SAVI -the index developed to compensate for the reflectance from soil (Huete, 1988), and WDRVI which was developed to increase the linearity with biophysical parameters (Gitelson, 2004 Table 4). Full-size DOI: 10.7717/peerj.5613/ fig-5 best proxy of GEP d (Table 4, Fig. 5). NDVI explained even 91% of GEP d variability of potato being the best predictor of GEP d . Inclusion of PAR into VI -based models, although successfully implemented in many studies (e.g., Gitelson et al., 2006;Wu et al., 2009;Rossini et al., 2010), did not improve goodness-of-fit of the linear regressions for any of the crops. Similar results have been reported also by Rossini et al. (2012) andSakowska et al. (2014) for alpine grassland ecosystems. Sakowska et al. (2014) hypothesized that this might be the result of a different response of plant photosynthesis to direct and diffuse radiation, but this cannot support our results, because all measurements were taken under sunny days. As stated by Rossini et al. (2012), although models including PAR and VI take into account variations related to changing incident radiation, they may not improve estimation of GEP due to higher light use efficiency (LUE) of plants at lower values of incident PAR and lower LUE at higher PAR (due to higher photoinhibition). Rossini et al. (2012) argued that higher photosynthetic efficiency at low PAR may be the result of two processes: (1) more diffuse light is penetrating more deep into the canopy, and (2) less photoinhibition on the top of the canopy which may reduce tendency towards saturation (Chen, Shen & Kato, 2009). It is well known that exposure of photosynthetic machinery of plants to strong light may result in inhibition of the photosystem II (PSII) activity, due to toxic effect of reactive oxygen species (Murata et al., 2007). Although this effect may be overcome by rapid and efficient repair of PSII (Aro, Virgin & Anderson, 1993), environmental factors, such as e.g., heat stress, which often occur during growing season, may inhibit the reparation of PSII and hence reduce efficiency of CO 2 uptake (Murata et al., 2007). This is probably why we observed a very weak correlation between GEP d and average daily values of incident PAR (Fig. 6). Although GEP d was increasing with increasing PAR during the vegetative phase of crop development, LUE became stable, while LAI and biomass of plants have been increasing continuously (Figs. 2 and 6). During the reproductive phase, both GEP d and LUE decreased with decreasing values of PAR d , but this effect is clearly related to progressive degradation of photosynthetic apparatus towards the senescence. Considering the above, we may hypothesize that the increasing GEP d observed during the vegetative phase of plant development cycle is mostly related to increasing biomass of plants and amount of photosynthetic apparatus rather than to increasing PAR. Another point is that the chamber measurements of CO 2 fluxes were taken on sunny days, when PAR was not a limiting factor for GEP. This is probably why we found a week correlation between GEP d and average daily PAR, which can indicate that over such a long seasonal time scales PAR is not the most relevant determinant of GEP d , although in shorter time-scales it may be more important (as indicated also by Sims et al., 2006).
One of the possible explanation, why the simple greenness indices considered in this study were moderately correlated with GEP d of winter crops might be that we considered all the 3-years period dataset together for all the crops, while it is clear that climatological conditions for all these years were extremely different. Especially, 2011-a year with a very dry spring, has disturbed the relationships between GEP d and VIs (mainly for winter crops) and although not shown, the same VIs-based models obtained specifically for 2012 and 2013 resulted in a much higher accuracy. However, our intention was to test the multiyear relationships between VIs and GEP d in order to develop more general models which can be applied to years characterized by different climatological conditions, still maintaining a reliable performance.
Another reason for higher uncertainties of GEP d estimations with VIs, observed specifically for winter crops, might be related to seasonal changes in canopy structure and biochemical traits which may modify the spectral response of plants over the growing season at different phenological phases of plant development (Asrar et al., 1984;Asrar, Kanemasu & Yoshida, 1985;Peng et al., 2017). The seasonal courses of GEP d and NDVI are nearly the same and are overlapping for potato (Fig. 7) and that is why the greenness indices (NDVI, SAVI and WDRVI ) are good proxies of GEP d for this crop (R 2 for GEP d vs. VIs relationships of 0.85, 0.90 and 0.91 for WDRVI, SAVI and NDVI respectively, see Table 4). However, in case of winter crops and spring barley these GEP d and NDVI seasonal courses did not correspond as well as they did for potato. We observed clearly a shift in seasonal courses of NDVI and GEP d for spring barley, winter rye and winter wheat, and peaks of NDVI occurred earlier than peaks of GEP d for all these crops (Fig. 7). It is well known that NDVI saturates at moderate-to-high biomass conditions (Gitelson et al., 2003), but this effect seems to be crop specific and may depend on the structure of the crop canopy. In spring barley NDVI was increasing together with LAI and GEP d until DOY140 (Fig. 7). After this day NDVI saturated and stabilized at around 0.8 until DOY165, although GEP d and LAI were already decreasing. For winter wheat, NDVI also followed LAI development during the vegetative phase (untill DOY140), but until DOY120 the increase rates of NDVI were slower and its value did not exceed 0.4 since the beginning of the growing season, while GEP d increased much faster during this period. Hence, we can hypothesize that GEP d vs. NDVI relationships might be different in the period between DOY80 to DOY120 and between DOY120 to DOY140, although we cannot confirm this due to not sufficient amount of data. After DOY140, NDVI saturated and stabilized again at around 0.8 until DOY165 and just after it begun to decrease slowly and the rate of decreasing was the same as the rate of GEP d decrease. Hence, in the reproductive/senescence phase of wheat development from DOY165 till the harvest, again GEP d vs. NDVI relationships were significantly (p < 0.05) different (Fig. 8). Much more complex analysis are related to winter rye, where NDVI reached maximum values of around 0.9 very soon after the beginning of the analyzed period and was quite stable until DOY140, although LAI and GEP d were continuously increasing until DOY145. When analyzing the winter crop data, specifically winter rye, one should take into account that this crop is sawn in the late September under the climatic conditions of the Central Europe and during warm winters it can continue to grow. In early spring, after the beginning of the growing period, the canopy of the crop may get very dense and green if it is growing under non N-limited conditions. In our case, the winter of 2012 was warmer than in 2011 and 2013, with a very warm December and January 2012 (Fig. 1), hence crops continued to grow during this time. Considering above, this may explain why NDVI saturated already in the late March/beginning of April 2012 at winter rye fields (Fig. 7). In 2011 and 2013 this effect was also observed, but beginning of the period when NDVI started to saturate occurred one month later (in April), due to longer and colder winters. For the first part of the growing season of 2012 for the winter rye, NDVI was not sensitive to changes neither in biomass, nor in GEP d (the same was observed in all the three analyzed years, data shown in (File S7). Between DOY160 and the harvest, changes of NDVI followed changes of GEP d , but again the GEP d vs. NDVI relationships were significantly (p < 0.05) different than those found for vegetative phase (Fig. 8). Similar kind of hysteresis were found in relationships between GEP d and LAI and GEP d vs. chlorophyll content in maize by Gitelson et al. (2014), as well as between reflectance in red and blue bands and greenness indices (e.g., NDVI ) vs. chlorophyll content in maize and soybean by Peng et al. (2017). The relationships between these variables were significantly different for vegetative and reproductive phases of maize and soybean development. In our study we found similar kind of differences between GEP d and NDVI in both analyzed winter crops, but not in spring barley and potato. Each of the crops has different height, LAI, canopy architecture, and different contribution of soil to canopy reflectance. As already indicated e.g., by Gausman et al. (1971) and Gausman, Rodriguez & Richardson (1976), and discussed by Peng et al. (2017), optical properties of leaf and canopies are greatly impacted by leaf structure and canopy architecture and that is why relationships between biophysical properties of crops and VIs are often different at different phenological phases of plant development (Asrar et al., 1984;Asrar, Kanemasu & Yoshida, 1985;Gitelson et al., 2014). In winter crops the canopy is very green, dense and more ''closed'' at the beginning of the growing season and hence reflectance from soil is minimized, while absorption in the red part of the spectrum is high. That is why NDVI saturates early in the season and it does not change much over the vegetative phase of crop development. Whereas in the reproductive/senescence phases, when the leaf structure and canopy architecture change and chlorophyll content is reduced, more light can penetrate deeper into the canopy and that is why the soil reflectance contribution is also higher. Due to differences in canopy architecture and leaf structure and different patterns of light absorption and reflectance by crops as well as different contribution of soil into overall reflectance of vegetation canopies at different phases of phenological development of plants, the relationships between GEP d and VIs of winter crops are different for vegetative and reproductive phases, as indicated also by (Peng et al., 2017;Gitelson et al., 2014). In potato, canopy is more open since the beginning of the growing season, and soil reflectance contributes significantly to the overall canopy reflectance. With the development of the green biomass, absorption of the red part of the spectrum is increasing, while soil contribution is decreasing. After the peak of biomass, during the reproductive and senescence phases, when chlorophyll pigments are degrading while leaves are folding, more light penetrates deeper to the canopy and again contribution of the soil reflectance to the overall canopy reflectance increases, whereas absorption in the red part of the spectrum decreases. Probably this effect can cause lack of hysteresis in GEP d vs. NDVI relationships for this crop, as the slopes of curves for this relationships in both phases of potato development are the same (Fig. 8).
In winter crops, LAI started to increase since the beginning of the growing season throughout the vegetative phase and GEP d followed changes in the crop biomass, although NDVI was already not sensitive enough to track changes in photosynthesis. The same LAI development was found in spring barley and potato, however NDVI followed changes in LAI and as discussed above, it was sensitive enough to track changes in GEP d of spring crops. After the peak of biomass, LAI of all the crops decreased slightly and stabilized at around 2-3 m 2 m −2 until the harvest (note, we are not investigating greenLAI, but total LAI ). This led to the conclusion, that LAI can be successfully used to overcome problems with greenness indices, which seem to saturate very early in the growing season in winter crops. By multiplication of NDVI -based VIs by LAI this specific issue of seasonal GEP d vs. VIs relationships can be overcome and model uncertainties can be reduced (see Table 4). The curves presenting the seasonal variations of the NDVI*LAI product are more closely related to seasonal changes in GEP d (Fig. 7) and this effect is not only restricted to winter crops, but can also be observed in case of spring crops. For both winter crops the slopes of the GEP d vs. NDVI relationships are significantly different (p < 0.05) between vegetative and reproductive phases and between the species, which can limit their application for accurate estimation of GEP d of these canopies over the entire season (Fig. 8). However, after multiplying NDVI by LAI slopes for the same relationships are much closer to each other (Fig. 8), and uncertainties of GEP d estimations are lower, as indicated in Table 4. Moreover, this approach does not change significantly slopes of the GEP d vs. NDVI relationships for spring crops (Fig. 8), although uncertainties of GEP d estimations for potato can be higher (Table 4).
However, in the more general, crop-combined dataset representing both winter and spring crops, the multiplication of NDVI and LAI led to an improvement of GEP d estimations (Table 4). NRMSE for the model fed with NDVI*LAI was about 14%, and it explained around 74% of GEP d variability independently from the crop species. The obtained accuracy is in the range which can promote this approach in remote sensing studies in order to overcome the hysteresis of GEP d vs. VIs relationship between vegetative and reproductive phases, which indeed may limit their applicability to predict photosynthesis of different crops over the entire growing season.
The limitation of this approach is that VIs are the remote sensing source of information which can be obtained from space, while space born products of LAI are result of (1) statistical models which quantify the relationships between LAI and canopy reflectance or VIs (Baret & Guyot, 1991;Verrelst et al., 2015;Linker & Gitelson, 2017), or (2) different radiative transfer models (Baret et al., 2007;Richter et al., 2012). These relationships between LAI and VIs are often canopy structure-and land-cover depended and are highly impacted by leaf angle distribution, vegetation clumping, optical properties of leaf and canopies (Goward & Huemmrich, 1992). What is more, different canopies may exhibit large variations in reflectance properties which can result in different values of VIs for similar values of LAI and other biophysical parameters (Pinty, Leprieur & Verstraete, 1993). There are few satellite based products of LAI (Zheng & Moskal, 2009) which are based not only on single vegetation indices such as EVI (Huete et al., 2002) or Reduced Simple Ratio (RSR, Brown et al., 2000), but also on linear or non-liner models including many vegetation indices which are used to estimate and map LAI at the landscape and global levels with Landsat satellites (Cohen, Maierspergerr & Tumer, 2003). Hence, we believe that this kind of data (LAI retrieved from space borne data) can also be used to mitigate the hysteresis effect described above and increase the accuracy of GEP d estimations for croplands. We can even speculate, considering that greenLAI and total LAI analyzed in this study are the same in the vegetative phase of plant development, while greenLAI is decreasing till ''zero'' in the senescence phase, although the total LAI remains relatively constant during this phase, that the product of NDVI and greenLAI retrieved from satellite data will improve estimations of GEP d of winter crops even better than total LAI used in this study. GreenLAI multiplied by greenness related VIs shall improve model accuracy both in the vegetative and reproductive phases of plant development cycle, due to reasons described above. Hence, we hypothesize that the overall hysteresis effect observed in relationships between greenness related VIs and GEP will be reduced. However, this effect will need to be studied in the future in more detail with the application of both greenLAI estimated at the ground and based on satellites data.
The proposed approach and empirical models developed in this study can be tested in other regions and for other C3 crops in order to verify the validity of our assumptions. Considering that our models were developed based on four different crops and three years characterized by different climatic conditions we believe that application of the proposed approaches and formulas can result in reliable estimation of daily GEP values of crops also for other regions with similar climate and crop management systems.

CONCLUSIONS
The analyzed multiyear relationships between GEP d and VIs showed that only in the case of spring crops GEP d can be estimated with a high accuracy and with an error smaller than 12% based on simple greenness indices (NDVI, SAVI, WDRVI ). The same kind of analyzes conducted for winter crops are much less accurate, and the error of GEP d estimation is higher than 18%.
The reason for the weaker correlation between daily GEP and VIs of winter crops may be related to hysteresis effect of this relationship found between the vegetative and reproductive phases of plant development cycle. We found that multiplication of greenness indices by LAI (which is much more sensitive to changes in biomass and GEP d of winter crops, specifically during the vegetative phase of their development) can mitigate such effect. The product of multiplication of LAI and VI has in most cases the same seasonality as GEP d and that is why it represents well the seasonal changes of gross CO 2 fluxes of croplands.
In order to propose as universal model as possible, we investigated the relationships between VIs and GEP d for the cereals-and crop-combined datasets, where both winter and spring crops, as well as cereals and potato were included. We found that there is no difference in GEP d estimations between these two kind of approaches and this general model based on approach where NDVI is multiplied by LAI, can be successfully applied for both winter and spring crops as well as for cereals and potato, while the error of this estimation is not higher than 14%. However, this approach underestimated GEP d of winter crops and potatoes, and overestimated GEP d of spring barley, but the rate of over-or under-estimations was not higher than 25%.