You could try local Moran's I index  for each year to identify the best search radius to use in Kriging intepolation. Read this article.

Article Quantitative monitoring of Aedes albopictus in Emilia-Romagn...

I would use a cross validation study to learn how prediction quality is related (if at all) to search radius.  The other aspect may be the geometry of the search neighborhood--should it be circular, elliptical rectangular?  You may also want to experiment with using a fixed number of points--such as the k-nearest neighbors and iterate on k. 

This is an empirical issue with little or no real theory available to make the selection.

Lastly, I would also use cross validation to compare your kriged surfaces with those from natural neighbor interpolation, as frequently the marginal improvement is small, at the expense of much effort to develop the kriging model.

Off the top of my head, I think it's better to use the same search radius for all the four years, using the average sampling distance.

I assume that the four "practical ranges" you have given refer to the variogram model, this is not the same thing as the choice of the search radius for kriging. In general when the distance from a data location to the interpolation point exceeds the variogram range that data location will have little influence on the interpolated value hence generally speaking you only want the search radius to be as large as or even less than the range of the variogram.

For example if all of the data locations in the search neighborhood are farther away from the interpolation point than the range of the variogram then the interpolated value is just the average of the data values

You need to look at a couple of other factors to determine an optimal (although that term is misleading) choice of the search radius. (1) How densely are the data locations placed around each interpolation point? In general if you use more than  the nearest 25 data locations then the others will have very small weights. (2) Will there be enough data locations inside the search neighborhood for a given search radius? If the data interpolation point is near the edge of the total set of data locations then you can end up with nearly all the data locations in the search neigborhood on just one side. The kriging variance will be  larger

One of the suggestions to use cross validation is a good one, depending on which software you use you should get multiple cross validation statistics (mean error, mean normalized squared error, fraction of errors greater than 2.5 (in absolute value).  A plot of the cross-validation errors (both normalized and non-normalized) is useful. You want the mean error close to zero, the mean normalized squared error close to one and the fraction of errors less than .05 (Chebyshev's inequality). The mean error can be misleading, if the the magnitudes of the data values are on the order of 1000 then an mean error of 5 is not large but if the magnitudes are on the order of 5 then a mean error of 1 is large. "Normalized" means divide by the kriging std deviation.   If you get any negative kriging weights then the search radius is too large

The theoretical assumptions for Moran's I are quite different from those for the kriging equations so that statistic is likely not relevant. If your software allows you to set a maximum number of data locations in the search neighborhood try using smaller numbers, 15 or less. If there are any interpolation points where there is no interpolation value it means that there were not sufficient data locations in the search neighborhood, i.e. when using cross validation. Specifying the range and/or the search radius to one decimal point is probably mis-leading, it won't make that much difference

If the fitness of your variogram model is good (meaning the model describes correctly the spatial behavior) than, in theory, the search radius should not be of much concern since the kriging weights will reflect your intended purpose.

Beware that limiting the search radius using kriging procedures may lead to image artifacts (usually wide concentric lines that seem to separate two areas in a way that defies common-sense).

On the event of the variogram model being bad than I would advise priorities. For instance if you have hard-data with distances in between them below 25 meters than it makes little sense to adjust a variogram model to predict distance far above. It would be better spent accurately predicting the first 25 meters (since the great majority of nodes will be within that distance) even if distances above will be bad.

Personally, if hard-data locations remain the same in all years and the anisotropy ratios are similar, I would use the same parametrization (search-radius, number of nodes, etc.) for all of them. The reason is that you intend to study the stability of your maps in four different years (which I assume that would be if they look alike or not) and I would rather not cause more uncertainty by playing with the parameters (it would be difficult to decide if the difference comes from the data or your interpretation of the data). This is, nevertheless, an interpretation and not a "by the manual" procedure.

Addendum

Given that the apparent variogram ranges vary so much over time have you thought about what might cause this, e.g. different weather patterns or something else?

Note that Inverse Distance Weighting is not a "derived" interpolator, i.e. it is not based on any statistical assumptions in contrast to kriging. IDW is also not an "exact" interpolator, i.e. you can't use it to obtain an interpolated value at a data location by leaving out the data value as you can do in kriging.

How do the variogram ranges compare with the distance across the region(s) of interest, much smaller, comparable or possibly even nearly the same? How did you determine the variogram ranges, i.e. how reliable are the values

We have used LASSO and cross-validation in our recent paper, which may be helpful. The paper is free online.

Zhang K. et al., 2014. Characterizing spatial patterns of airborne coarse particulate (PM10-2.5) mass and chemical components in three cities: The Multi-Ethnic Study of Atherosclerosis, Environ Health Perspect 122(8): 823-830.

Dear Moshood,

As per the project, it's better to use the same search radius for all the four years, using the sampling points, then only, it is useful for map overlay and comparison, and can be found changes.

Best regatrds.

M.Palaniyandi

kriging is called BLUE “the best linear unbiased estimator”

Your question seems to imply that the spatial correlation is temporally dependent but you didn't compute or use a spatial temporal variogram model.

Are the data locations the same for each of the four years? If not, was there a good reason they are different? If they are the same then the search neighborhood parameters should be the same.

You say that "the spatial pattern is best described by an exponential variogram", that is nonsense and meaningless.

You give the practical ranges with one decimal point precison, this is nonsense. It is not possible.

As I have already pointed out, you are confusing the range of the variogram with the radius of the search neighborhood.

There is no justification for claiming "unbiasedness" if you use the same search radii for the searches.

What you need to do is attempt to discern how sensitive the kriging results are to changing the radii, i.e. use cross validation. Note that there are several different statistics that might be used to quantify the sensitivity but if you attempt to use them to optimize the choice the different statistics may conflict.

Azouzi's comment is meaningless because it does not recognize that there are several different forms of kriging, i.e with slightly different statistical assumptions. Moreover the grammar is wrong, you would have to say that the kriging estimator is a BLUE

Alessandro's suggestion is wrong because the statistical assumptions pertaining to variograms are completely different from those for Moran's I index. Moran's I is empirical whereas a variogram is a function (conditionally negative definite).