## 2 Materials and Methodology
### Data Analysis
#### 2.3.2 Trends and Variability in GPP  
Long-term changes in GPP were assessed for the period of 1982-2016. The magnitude of change was determined using [PERSON]'s slope estimator, while the statistical significance of these changes was assessed using the non-parametric Mann-Kendall (M-K) statistics. The M-K test does not assume a priori distribution of the data and is less sensitive to the presence of outliers and anomalous values in the data [69]. For this reason, it has been widely applied in different environmental applications [70; 71; 72; 73]. The M-K test's S statistic [74; 75] was computed as:  
\[M-K\left(S\right)=\sum_{k=1}^{n-1}\sum_{j=k+1}^{n}sgn\left(x_{j}-x_{k}\right) :\left(j>k\right), \tag{1}\]  
where \(x_{j}\) and \(x_{k}\) are the GPP values in years \(j\) and \(k\)\(\left(j>k\right)\), \(n\) is the time series length, and \(sgn(\theta)\) is the sign function expressed as:  
\[gn\left(x_{j}-x_{k}\right)=\begin{cases}+1\big{(}if\big{(}x_{j}-x_{k}\big{)}> 0\big{)}\\ 0\big{(}if\big{(}x_{j}-x_{k}\big{)}=0\big{)}\\ -1\big{(}if\big{(}x_{j}-x_{k}\big{)}<0\big{)}\end{cases}, \tag{2}\]  
When \(n\geq 8\), the \(S\) statistic has a normal distribution, and its average and variance can be computed as:  
\[E(S)=0, \tag{3}\]\[\text{Var}(S)=\frac{n\ (n-1)(2n+5)-\sum_{i=1}^{m}t_{i}(t_{i}-1)(2t_{i}+5)}{18}\, \tag{4}\]  
Equation (4), which was used to calculate the variance in M-K (S), was based on the idea that there are tied values in the series (i.e., the same consecutive data) for modifying the variance, where m denotes the number of tied groups and is the number of data values in the \(i\)-th group. The \(Z\) statistic of the M-K test can be computed as:  
\[Z=\left\{\begin{array}{c}\frac{s-1}{\sqrt{\text{var}(s)}}(s>0)\\ 0(s=0)\\ \frac{s+1}{\sqrt{\text{var}(s)}}(s<0)\end{array}\right., \tag{5}\]  
where the \(Z\) statistic was applied to examine the existence or non-existence of statistically significant trends at \(\alpha\) significance level (a positive \(Z\) suggested a rising trend in the time series and a negative \(Z\) denoted a declining trend in the time series). Herein, the null and alternative hypotheses were tested at a significance level of 95% (\(p<0.05\)).  
To quantify the change ratio of the GPP trend (i.e., the change per unit time), the non-parametric [PERSON]'s slope estimator was applied [76; 77]. The merit of using this method is that it is not sensitive to the presence of extreme values; it can be computed as:  
\[median(bi)=\frac{x_{j}-x_{k}}{j-k},where\ j>k, \tag{6}\]  
where \(x_{j}\) and \(x_{k}\) are the GPP data values for years \(j\) and \(k\) (\(j>k\)). The Sen's slope \(bi\) is the median of \(n\) values of \(b\). The \(n\)-values of \(b\) are ranked from the smallest value to the largest, and the Sen's slope is defined as:  
\[Q=\left\{\begin{array}{c}b(N+1)/2\ \  
where ERI is the ecosystem's GPP resilience index, \(GPP_{i}^{o}\) is the observed GPP in the driest year, and \(GPP_{i}^{(\tau)}\) is the value of the de-trended value of GPP in the same year. The values of the ERI could be classified into four categories, which were predetermined for crop resilience as follows: slightly non-resilient (0.9 < ERI < 1), moderately non-resilient (0.8 < ERI < 0.9), and severely non-resilient (ERI < 0.8) [80,82]. The driest year for each pixel was defined using the lowest annual SEDI over the study period.  
Finally, the Pearson correlation (\(r\)) was used to examine the relationships between the SEDI and GPP as well as to detect the consistency among the GPP models (i.e., the spatial cross-correlation analysis) in terms of trends and their resilience to severe drought.