Are Kriging's residuals (i.e. $Z-\hat{Z}$) spatially independent? Suppose data ${Z(s_i ):i=1, ..., n}$ are observed at spatial locations  ${s_i :i=1, ..., n}$. To carry out the spatial prediction (predict un unknown $Z(s_0)$ at a known location $s_0$) we can use a method like ordinary Kriging.
I've used cross-validation for Kriging to create residuals to compare actual versus model’s predicted values. I know that for these residuals, RMSE (root mean square error) should be as small as possible. But my question is that is there any spatial dependence in these residuals?
Should Kriging's (or any other spatial prediction method's) residuals be independent? If the answer is YES, can I use this to compare two different spatial prediction methods? (I mean, the more the method remove the spatial dependency in residuals, the better it is).  
 A: $$
\newcommand{\E}{\mathbb{E}}
\DeclareMathOperator{\cov}{Cov}
\newcommand{\zhat}{\widehat{Z}}
$$
I'll try and answer the title of your question, about spatial independence of the residuals.
For ordinary kriging, the covariance matrix of $Z = f + \epsilon$, where $f$ is the spatial (Gaussian) process and $\epsilon$ is iid noise, decomposes into the spatial covariance $K$ and a diagonal term $\sigma_n^2 I$. We can write $\cov(Z) = \Sigma = K + \sigma_n^2I$.
The predicted values are obtained as the posterior expectation of the underlying process
\begin{equation}
\begin{split}
\zhat &= \E(f \mid Z) = \cov(f, Z) \cov(Z)^{-1} Z \\
&= K \Sigma^{-1} Z\,.
\end{split}
\end{equation}
The prior covariance of the residuals can be derived analytically:
\begin{equation}
\begin{split}
\cov(Z - \zhat) 
    &= \cov(Z) -\cov(Z, \zhat) -\cov(\zhat, Z) + \cov(\zhat) \\
    &= \cov(Z) -\cov(Z, K \Sigma^{-1} Z) -\cov(K \Sigma^{-1} Z, Z) + \cov(K \Sigma^{-1} Z) \\
    &= \Sigma -2 K \Sigma^{-1} \Sigma + K \Sigma^{-1} \Sigma \Sigma^{-1} K \\
    &= \Sigma - 2K + K\Sigma^{-1}K \\
    &= \sigma_n^2I - K + K\Sigma^{-1}K
\end{split}
\end{equation}
which is not going to be diagonal.
Another way to derive the covariance of the residuals is to use the fact that for any random variables $X$ and $Y$, $Y - \E(Y \mid X)$ and $h(X)$ are uncorrelated from all measureable functions $h(X)$ of $X$.
Choosing $Y=f$, $X=Z$, and $h(Z)=Z-\zhat$, we get that
\begin{equation}
    \cov(f - \zhat, Z - \zhat) = 0\,.
\end{equation}
So
\begin{equation}
\begin{split}
    \cov(Z-\zhat, Z-\zhat) 
        &= \cov(\epsilon + f - \zhat, Z-\zhat) \\
        &= \cov(\epsilon, Z-\zhat) + \underbrace{\cov(f - \zhat, Z-\zhat)}_0 \\
        &= \cov(\epsilon, \epsilon+f-\zhat) \\
        &= \underbrace{\cov(\epsilon)}_{\sigma_n^2 I} + \underbrace{\cov(\epsilon,f)}_{0} - \cov(\epsilon, \zhat) \\
        &= \sigma_n^2 I - \cov(\epsilon, K \Sigma^{-1} (f-\epsilon)) \\
        &= \sigma_n^2 I - \sigma_n^2 K \Sigma^{-1}
\end{split}
\end{equation}
which again does not look like it will be diagonal.
These two results can be shown to be the same, i.e. $K - K\Sigma^{-1}K = \sigma_n^2 K \Sigma^{-1}$, by left-multiplying this identity by $K^{-1}$ and then $\Sigma$.
