The biggest issue I have with the paper is trying to figure what was done in the methods.   Obviously, we need this kind of comparison in climatology and it is extremely important.   I didn't have any issues with the scietific prose in the paper, but I would like the authors to be much more pedagogical in their methodological exposition so that folks can discern what has been done.   I did not get anything out of the brief descripton of the four methods in Section 2.   Hence, I was hoping that Section 3 would alleviate these concerns; alas, it did not.    

In Section 3, some things are not stated that need to be:   are the X and Y processes independent (I think so)?    It seems you are assuming a constant mean in time t (which is not likely true) and that the covariance function of the spatial fields at each time have the same structure (I can buy this).  I also want to know if you are assuming that the fields are Gaussian.

In testing for whether the means of the two processes are the same, why would we not just look at the average (over all spatial locations and times) and use asymptotic normality to test whether these X minus Y averages have a zero mean?    This just works with differences \Delta(t,s)=X(t,s)-Y(t,s).  Then you don't have to assume the mean is constant....it subtracts to zero under the null.  You can easily estimate the variances of the average \Delta value assuming a null that the two fields have the same covariance structure.     This seems to be the fundamental way to handle the two sample equality issue in general abstract spaces.    I'm guessing that what you've done can be justified, but it would seem that I have to go to your past papers to dig this up.   I just have this uneasy feeling that the EOF approach is needlessly complicated.  I also can't rationalize why I need to use the data from times 1 to k in various places.   Seems I should use all N times once.   

You've also got some key typos and omissions here:   V_\psi(L) versus V_\psi is confusing and you've never told me the distribution (asymptotic or not) of TS1.    I presume this is chi-squared, but I was left wondering (Ditto TS2).    And is V_\psi some sort of covariance estimate?   If so, please note.

It would seem to me that we want to test whether the means are the same and the covariances are the same in tandem.   Not either the mean is the same or the covariance is the same separately, but to test both in tandem.   So why not set \Delta(t,s)=X(t,s)-Y(t,s) and work with these differences as above.   If the means are the same, the mean of the \Delta process is identically zero at all times and spatial locations.  Then we could stack the \Delta(t,s) in a giant vector --- call it V --- over all spatial locations and time points.   Now if we could get the covariance matrix of all components in this giant vector --- call it \Sigma --- we would just look at \Sigma^{-1/2}V.   This quantity would be composed of IID N(0,1) variates if the original fields are Gaussian.   And it is easy to test whether data is IID N(0,1) by a plethora of methods (QQ plots, Kolmogorov-Smirnov, chi-squared tests, etc).   To estimate \Sigma, grab your favorite space-time covariance estimators to estimate both the X covariance and the Y covariance structures in time and space.   Call these estimates \Sigma_X and \Sigma_Y, respectively.  Let \Sigma^*=(\Sigma_X+\Sigma_Y)/2 be the common estimate under the null that the two processes have the same covariances and are independent.   Now just use Cov(\Delta(t,s), \Delta(t^\prime,s^\prime))= 2 times the corresponding entry in the matrix \Sigma^*.    Then I think it's game over:  you've tested both hypotheses at the same time.  

I can spell this out in more detail if you need it.