## 1 Introduction  
Satellite observations of vegetation greenness have been used as a means to characterize amount, rate, direction, location, timing, drivers and consequences of changes in ecosystem structure and function at various spatio-temporal scales [1-6]. Seasonal and interannual vegetation dynamics and phenological patterns at the ecosystem level (e.g., timing and rate of green-up and senescence of vegetation classes, and amplitude and duration of growing season) constitute one of the key driving variables for modeling and monitoring of terrestrial ecosystems [7-9]. In response to natural and/or anthropogenic effects, the scientific community is increasingly interested in deriving information about and tracking changes in ecosystems from remotely sensed data, in order to better understand the natural and human-induced processes of changes in land cover (LC) and land use (LU) over time and space, their implications for biogeochemical cycles, and ways of interventions through management, planning, and policy to ensure sustainability of ecosystem goods and services [10-12]. Moderate Resolution Imaging Spectroradiometer (MODIS) as an instrument on board NASA's Terra and Aqua platforms for remote sensing of the atmosphere, oceans and land surfaces provides vegetation indices (VIs) more accurately than Advanced Very High Resolution Radiometer (AVHRR) as a highly evolved successor to AVHRR [13]. The two MODIS VIs of Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) on a scale of minus one (-1) to plus one (+1) are spectral measures of the amount, relative greenness, phenological characteristics, and biological productivity of observed vegetation present on the ground, with a global coverage every one to two days at three spatial resolutions of 250 m, 500 m, and 1 km [14, 15]. NDVI and EVI are calculated using the following equations [13]:  
\[\mathrm{NDVI}=\frac{NIR-R}{NIR+R}\,, \tag{1}\]  
where _NIR_ and \(R\) stand for the spectral reflectance measurements acquired in the red and near-infrared regions, respectively. \[\mathrm{EVI}=G\,\frac{NIR-R}{NIR+C_{1}R-C_{2}B+L}\,, \tag{2}\]  
where _NIR_, \(R\), and \(B\) are reflectances in the near infrared, red, and blue bands respectively; _C_\({}_{1}\) and _C_\({}_{2}\) are aerosol resistance coefficients; \(G\) is a gain factor, and \(L\) is the canopy background adjustment that addresses nonlinear, differential NIR and red radiant transfer through a canopy. The coefficients adopted in the MODIS-EVI algorithm are \(L\) = 1, _C_\({}_{1}\) = 6, _C_\({}_{2}\) = 7.5, and \(G\) = 2.5 [16]. One of the most common VI-based approaches to quantify net and gross primary productions over time is through the following equation in process-based models [17-20]:  
\[\mathrm{NPP}=\epsilon_{n}\times\mathrm{FAPAR}\times\mathrm{PAR}\;\;\mathrm{ or}\;\;\mathrm{GPP}=\epsilon_{g}\times\mathrm{FAPAR}\times\mathrm{PAR}\;\;, \tag{3}\]  
where  
\[\mathrm{FAPAR}=a+b\times\mathrm{NDVI}\,, \tag{4}\]where NPP and GPP refer to net and gross primary productivity, respectively; PAR is the incident photosynthetically active radiation (MJ m\({}^{-2}\)) for a given time period (e.g. hour, day or month); FAPAR is the fraction of PAR absorbed by vegetation canopy and is considered to be linearly related to NDVI; and \(\varepsilon_{n}\) and \(\varepsilon_{g}\) are the light use efficiencies (g C MJ\({}^{-1}\) PAR) in the calculation of NPP and GPP, respectively. \(\varepsilon_{n}\) is the product of the potential light use efficiency (\(\varepsilon_{0}\)) and reduction factors scaled to 0 to 1 that reflect environmentally limiting and reducing conditions. \(a\) and \(b\) are coefficients in the linear regression model. Errors in VI time series are often caused by atmospheric and ground conditions (e.g. cloud and snow cover), and sensor problems (e.g. sensor drift, and changes in sensor view angle), thus creating irregularly low VI values or data gaps in time series.