3.1. Exploratory Data Analyses
According to NPP time series (Fig. 2) for different regions, it can be seen that the NPP has monthly fluctuation. In addition, the time series looks stationary during the data record as they do not show a significant increasing or decreasing trend.
The NDVI times series are also shown in Fig. 3 for different bioclimatic regions. Almost all time series have no significant change, except for the Khazari (in Jan 2008) and Semi-Steppe (in Jan 2008) regions which indicate a breakpoint in the time series structure.
To show the seasonal characteristics of NPP and NDVI, boxes plot of monthly NPP and NDVI shown in Figs. 4 as well as 5, respectively. Each box plot contains annual data for each month during record period. We can observe in the NPP time series that the seasonal patterns in Semi-Desert, Arid Forest, Steppe, and Semi-Steppe regions are similar, while Baluchi and Khazari regions show different seasonal patterns which due to the different climate of these areas.
The box plot of the NDVI index shows that the seasonal patterns in the Arid forest, Steppe, and Semi-Steppe regions were similar. Again, the Baluchi and Khazari regions have different seasonal patterns. Semi-Desert region did not have distinct seasonal patterns (Fig. 5). The seasonal average NPP and NDVI time series are also provided in Tables 1 as well as 2, respectively.
As it is shown the average NPP in Khazari region is highest in May. This amount is observed for the Baluchi region in February, for the semi-desert region in March, for the steppe region in April, for the semi-steppe region in February, and for the Arid forest region in April (Table 1). In cold regions, the maximum NPP is observed in spring, and in warm regions the maximum NPP is observed in winter which is due to the presence of suitable conditions for plant growth in these seasons.
The minimum amount of NPP is observed in Iran's Turani flora regions (semi-deserts, steppe, semi-steppes, and dry forests) in early summer as well as July that related to summers hydrated deficits duration, that leds to dry in farmlands vegetation biomasses (zoffuli et al., 2008)., in the Khazari region in winter and February related to vegetation inaction, and within Baluchi region in late spring and May due to dryness and high temperature in this month (Table 1).
Table 1
Monthly average NPP (kg/ha) during 2000 to 2016 for different bioclimatic regions (Bold values show maximum and italics are the lowest NPP (kg/ha) for each region)
| Jan | Feb | Mar | Apr | May | Jun | July | Aug | Sep | Oct | Nov | Dec |
Khazari | 45.43 | 51.98 | 70.81 | 106.98 | 135.2 | 126.09 | 102.74 | 76.80 | 95.50 | 84.33 | 61.65 | 45.57 |
Baluchi | 39.30 | 40.42 | 32.84 | 22.82 | 16.60 | 16.21 | 17.24 | 18.87 | 20.55 | 23.87 | 31.55 | 34.94 |
Semi-desert | 6.59 | 7.05 | 7.65 | 7.23 | 5.71 | 4.32 | 3.92 | 4.42 | 5.53 | 6.68 | 6.97 | 6.50 |
Steppe | 20.32 | 22.16 | 25.93 | 28.37 | 23.68 | 18.04 | 15.64 | 17.18 | 20.78 | 23.00 | 22.46 | 20.37 |
Semi-Steppe | 43.40 | 48.03 | 47.96 | 44.84 | 39.37 | 31.14 | 26.73 | 27.51 | 31.46 | 35.46 | 40.17 | 39.87 |
Arid forest | 39.63 | 44.80 | 59.45 | 76.20 | 64.18 | 34.25 | 22.35 | 24.93 | 37.34 | 43.04 | 44.30 | 39.42 |
The average NDVI from 2000 to 2018 for different months Table (2) showed that the Khazari regions, Semi-Steppe and Arid Forests have a distinct seasonal pattern. This attributed to the great variations in rainfall as well as heat among the seasons in these areas. In the semi-desert region, there is no difference in the amount of NDVI in different months, may be due to the limited moisture of plant growth throughout the year. The seasonal pattern in Baluchi and Steppe regions are not significant either because the rainfall in the wet season is very low and does not create a big difference between different seasons.
Table 2
Monthly average NDVI during 2000 to 2018 for different bioclimatic regions (Bold values show maximum and italics are the lowest NDVI)
| Jan | Feb | Mar | Apr | May | Jun | July | Aug | Sep | Oct | Nov | Dec |
Khazari | 0.25 | 0.25 | 0.31 | 0.41 | 0.46 | 0.45 | 0.45 | 0.43 | 0.43 | 0.43 | 0.35 | 0.30 |
Baluchi | 0.08 | 0.08 | 0.09 | 0.07 | 0.07 | 0.07 | 0.06 | 0.07 | 0.07 | 0.07 | 0.07 | 0.08 |
Semi-desert | 0.08 | 0.08 | 0.08 | 0.08 | 0.08 | 0.08 | 0.08 | 0.08 | 0.08 | 0.08 | 0.08 | 0.08 |
stepp | 0.10 | 0.10 | 0.11 | 0.12 | 0.10 | 0.10 | 0.10 | 0.10 | 0.10 | 0.10 | 0.10 | 0.10 |
Semi-Stepp | 0.14 | 0.16 | 0.18 | 0.18 | 0.16 | 0.14 | 0.13 | 0.13 | 0.13 | 0.13 | 0.13 | 0.14 |
Arid forest | 0.16 | 0.18 | 0.25 | 0.32 | 0.28 | 0.21 | 0.18 | 0.18 | 0.17 | 0.18 | 0.19 | 0.18 |
3.2. Time Series Modeling
3.2.1. Autocorrelation Structure
Examining the ACF and PACF functions makes it possible to understand how the time series behaves in stationarity and seasonality. If ACF decreases gradually, the time series is non-stationary; if ACF decreases suddenly, the times series is constant. Before developing times series models, the behavior of ACF of NPP and NDVI time series are examined to compare the different stationary and non-stationary characteristics of them.
For NPP times series, in lag-1autocorrelation, the autocorrelation value has high in all regions (0.6 to 0.8). That high lag-1autocorrelation for the Khazari and Semi-Steppe regions is most probably due to relatively continuous seasonal rainfall throughout the year and, despite the low annual rainfall (around 100 mm) in Baluchi, Semi desert and Steppe regions, they have high autocorrelation (about 0.7). In Baluchi region it is probably due to the presence of rich annual species and in Steppe region it is due to the existence of irrigated agricultural plains in some regions, such as Khuzestan province. In Semi desert region it can be due to the uniform environment, permanent lack of water, and vegetation adapted to these conditions. The yearly precipitation of Arid Forest area is high (around 600 mm) and it has suitable conditions for the growth of plants, but it has a lower autocorrelation value compared to drier regions which this can be due to large differences in temperature and precipitation between seasons (Fig. 6a, b, c, d, e, f). Considering that the autocorrelation value has suddenly decreased in all regions, it indicates the stability of NPP time series in all studied regions. All study regions have a 12-month periods regular seasonal patterns (Fig. 6a, b, c, d, e, f). In the semi-steppe and Bauchi regions, there is also a strong autocorrelation in 6-month periods, that is due to the growth conditions of plants due to suitable temperature and autumn rains (Fig. 6b, e).
The Arid Forest region has the different autocorrelation pattern. Many changes in this region can be due to climatic variations and the destruction of forest trees by humans and diseases such as Loranthus (Javanmiri pour et al. 2022) (Fig. 6f).
For the NDVI time series, in 1-month lag the ACF was 0.55 to 0.83. The Khazari regions has maximum autocorrelation (0.83) and Semi steppe region has minimum autocorrelation (about 0.55). In the Semi-Desert region, despite the inadequate environmental conditions for plant growth, ACF is higher (about 0.8) than the Semi-Steppe and Arid forests regions, which can be due to the presence of plant species adapted to dry environments (Fig. 7a, c, e, f). Since the ACF has suddenly decreased in Steppe, Semi steppe and Arid forest regions, it shows the stability of NDVI time series in this regions, and the ACF decreases gradually in Khazari, Baluchi and Semi desert was indicating non- stability (Fig. 7a, B, C, D, E, F). Time series of the Khazari, Baluchi, steppe, Semi-Steppe and Arid forests regions had periodic changes of 12 months. This is indicating that the annual cycle of vegetation phenology affects the NDVI. The semi-desert region had weak seasonal behavior which can be due to the very low vegetation cover in this area (Fig. 7a, b, c, d, e, f).
Comparing the lag times from 1 to 24 month for different regions shows that the NPP variable has a stronger seasonality than the NDVI (Fig. 8). In NPP, the maximum value of ACF is observed at the 12-month lag time while for NDVI it is observed at the 1-month lag. However, the ACF value at lag12 is also high for NDVI time series indicating significant seasonality.
Figure 8 also indicates that the range of autocorrelation coefficients is higher for NDVI time series comparing NPP time series which shows a distinct diversity between bioclimatic regions.
3.2.2. SARIMA Models
As the ACFs of NPP and NDVI show significant seasonality, the Seasonal ARIMA model seem to be appropriate. To fit the SARIMA model to the NPP and NDVI time series, regarding the seasonal non-stationary behavior of the time series, seasonal differentiation (D = 12) is considered. For the all-time series, different models with different magnitudes of non-seasonally and seasonally autoregressive as well as moving averages are investigated with and without non-seasonal differentiation. Finally, the best model was selected regarding minimum AIC criteria and no autocorrelation in the residuals the time independence of the residuals is checked using the ACF diagram in residuals. Normality of residual is also controlled through Q-Q plot diagram for NPP (Fig. 9) and NDVI (Fig. 10). Given model adequacy, the best SRIMA models for NPP and NDVI in each region is selected (Table 3).
Table 3
The selected SARIMA models for NPP and NDVI time series in all regions
| NDVI series | NPP series |
Khazari | SARIMA(0,1,1)(0,1,1) 12 | SARIMA(1,1,1)(1,1,1)12 |
Baluchi | SARIMA(1,1,1)(0,1,1) 12 | SARIMA (1,1,1)(0,1,1) 12 |
Semi-Desert | SARIMA(0,1,1)(1,1,1) 6 | SARIMA (1,1,1)(0,1,1) 12 |
Steppe | SARIMA(1,1,1)(0,1,1) 12 | SARIMA (1,1,1)(0,1,1) 12 |
Semi-Steppe | SARIMA(0,1,0)(0,1,1) 12 | SARIMA (1,1,1)(0,1,1) 12 |
Arid forest | SARIMA(0,1,0)(0,1,1) 12 | SARIMA (1,1,1)(0,1,1) 12 |
In order to check the efficiency of NPP and NDVI time series perdition, the error criteria are calculated (Table 4). These criteria show that the models for the NPP time series have more accuracy than NDVI time series based on RMSE, R2, MRE, and CE criteria, while based on the ME criteria, the models perform better for NDVI time series. Given that the difference between the minimum and maximum values in the variable NPP is greater and the ME criterion cannot weight the minimum and maximum values, this criterion is not suitable for the NPP time series models. The selected models for Baluchi and Steppe regions which have less climate variability, have higher accuracy than the more humid regions (Khazari, Semi-Steppe, and arid forest). For Semi desert region, the selected model for NPP time series have more accuracy than the NDVI time series (Table 4).
Table 4
Errors criteria for selected models in all region
| ME | RMSE | R2 | MRE | CE |
NPP | NDVI | NPP | NDVI | NPP | NDVI | NPP | NDVI | NPP | NDVI |
Khazari | 3.67 | 0.05 | 0.12 | 0.18 | 0.87 | 0.63 | 0.02 | 0.13 | 0.84 | 0.12 |
Baluchi | -1.23 | 0.01 | 0.13 | 0.23 | 0.91 | 0.99 | -0.10 | 0.04 | 0.85 | 0.15 |
Semi-Desert | -0.11 | 0.003 | 0.12 | 0.38 | 0.95 | 0.46 | -0.02 | 0.04 | 0.93 | -2.03 |
Steppe | -0.55 | 0.01 | 0.17 | 0.005 | 0.93 | 0.82 | -0.02 | -0.03 | 0.95 | 0.66 |
Semi-Steppe | 0.62 | 0.01 | 0.14 | 0.08 | 0.87 | 0.54 | -0.01 | 0.09 | 0.90 | 0.31 |
Arid forest | 0.61 | 0.04 | 0.12 | 0.22 | 0.89 | 0.46 | -0.04 | 0.18 | 0.84 | 0.15 |
3.2.3. Out of-sample Forecasting
Following specifying model as well as evaluating variables, next step is to evaluate the model's capability in out-of-sample forecasting between 2015 to 2016 for NPP and between for 2017 to 2018 for NDVI. For the NPP variable, the data from 2015 to 2016 were used to evaluate the model's forecasting. In the Khazari region, the model is able to forecast the changes both in trend and quantity well, and only in May 2015, there is a significant difference in the forecasted trend and value, which may be due to the indiscriminate harvesting of wood and other unpredicted environmental factors. The proposed model in the Balochi region also forecasts the changes well, and only in June 2016, both in terms of trend and value, there is a significant difference between observations and forecasting. In the semi-desert region, the model is able to forecast the changes well regarding trends and values, and all points are within the 95% significance range. In the steppe, Semi-steppe, and Arid forest regions, the models have forecasted the NPP values and trends very well, and all points are within the 95% confidence range. Generally, the selected NPP time series models can correctly forecast the NPP changes in all regions (Fig. 11).
We show the monthly observed and forecasted NDVI time series from January 2017 to December 2018 in Fig. 12. This figure shows which elected models to NDVI times series in the Khazari region performs very well in forecasting NDVI and its tendency, as forecasts are within the 95% confidence interval. In the Baluchi region, only in September 2017, the observed value was outside the 95% confidence interval of the model. The proposed models in Semi-Deserts, Steppes, Semi-Steppes, and dry forest regions have predicted the tendency and data values very well, and the observed data were completely within the range of 95% predicted by the model. Muti et al. (2019) also stated that SARIMA modeling performed better in predicting dry seasons, where the variability of climate is less, and predicts variable values in wet seasons, where annual changes are higher (Muti et al. 2019). Fernández Manso et al. )2011(also and stated that the prediction of the selected models for the NDVI data series using SARIMA was acceptable.
The performance criteria of the models in forecasting NPP and NDVI are presented in Table (5). The results of these criteria showed that for the NPP variable, the selected models for the semi-desert region performs better than other regions (R2 = 0.94, RMSE = 0.12, ME = 0.05, MRE = 0.01, and CE = 0.91) and other regions also have acceptable accuracy. For the NDVI index, the selected models for the Arid forest region performs the best (R2 = 0.87, RMSE = 0.22, ME = 0.044, MRE= -0.04, and CE = 0.86), and Khazari, Baluchi, and Steppe have acceptable accuracy, but the selected models for semi-desert and semi-steppe regions do not show high accuracy. In general, the results showed that time series modeling has a suitable capability for forecasting NDVI and NPP, and in all regions, the forecasted values are within 95% confidence.
Table 5
Errors criteria for forecast models in all regions
| ME | RMSE | R2 | MRE | CE |
NPP | NDVI | NPP | NDVI | NPP | NDVI | NPP | NDVI | NPP | NDVI |
Khazari | 3.67 | 0.05 | 0.13 | 0.19 | 0.86 | 0.62 | 0.02 | 0.14 | 0.85 | 0.11 |
Baluchi | -3.23 | 0.01 | 0.12 | 0.19 | 0.91 | 0.57 | -0.16 | 0.04 | 0.83 | 0.15 |
Semi-Desert | 0.05 | 0.003 | 0.12 | 0.38 | 0.94 | 0.43 | 0.01 | 0.04 | 0.91 | -2.03 |
Steppe | -0.67 | -0.01 | 0.17 | 0.004 | 0.92 | 0.80 | -0.03 | -0.03 | 0.93 | 0.60 |
Semi-Steppe | 0.62 | 0.013 | 0.14 | 0.02 | 0.85 | 0.38 | -0.01 | 0.08 | 0.87 | 0.24 |
Arid forest | 0.61 | 0.044 | 0.11 | 0.22 | 0.89 | 0.87 | -0.04 | 0.18 | 0.85 | 0.86 |