Sounds like a white noise (independent): in this case, it is indeed very difficult to interpolate (no spatial coherent structure).

Are you sure of your empirical estimation of the semivariogramm?

Maybe the best in this case, is to generate the same type of noise (with respect its histogram) to create value in-between ...

Please check using spatial autocorrelation. If it is the highest, it means this real case of invariance. Use either Type Qualifier or gain insight and find solutions from the following paper:

http://arxiv.org/pdf/1408.0210v1.pdf

Dear Stefan,

I suppose that you are facing a pure nugget effect. However, first of all is that with 33 observation points, having a cluster distribution, I think that you cannot construct a robust experimental variogram and then a robust model. In this condition the only way to provide a spatial model (i.e. interpolation) is to use deterministic simple algorithms for example a global and/or local trend surface analysis i.e., first or second order regression surface. Standard local and global regression interpolation algorithms are available in different packages, i.e. ArcGIS.

You should check the data, it may be that there has been an erroneous data entry, as if a Variogram is flat then it means that the trend is perfectly random and  this when making these studies  seems very unlikely

1. Since there is no variance or rather trend over distance, Toblers first law of geography (the nearer the more similar) is not fulfilled. The interpolation cannot be done based on a regression of values and distance (which is the variogram). Since there is no dependence of the values on the distance the mean value is the right result for all points in space when you assume that the data points has to map a trend. In this case spline or linear interpolation would produce a result with different values which are not correct.

2. The data set could consist of minimal information. Imagine a hilly landscape, were the data points are on the hilltop and in the valleys and the hills have the same altitude. You would have only two values which do not show a distance dependency (top, bottom, top, bottom, ...). In this case a spline would produce a perfect result, while kriging based on a semivariogram would fail. In this case spatial autocorrelation would be negative.

3. I agree with Giancarlos comment on a nugget effect. Only two points at the same location with very different values can destroy your model. You should check if you have wrong data.

If it is perfectly random the experimental variogram points should be more or less over the sill (assuming you're using data variance as sill).

The phenomena which causes an effect like that are uncommon but they exist (for instance vertical direction in a well log - rapid change of geological layer). In any case I agree with previous colleagues that said you should check your data for any error that may have occurred since that is not what I would expect from most phenomena.

Notice that using just a nugget effect as variogram model for geostatistics is the same that saying that all samples have the same weight no matter the distance to the point which is being interpolated/estimated.

Adding just one more suggestion to Giancarlos answer you can try to optimize a map by using stochastic procedures (still within geostatistics) such as sequential gaussian simulation or direct sequential simulation.

This is, for instance, common in procedures such as seismic inversion.

It could ne hwlpful to know a bit more on the nature of your data an on the spatial distribution of your samples.

Maybe you do not have samples close enough to capture the spatial variability of the process you want to study.

Moreover, perhaps you could stratify the samples beforehand.

33 data locations is not very many, the total number of paired differences is fixed by the number of data locations. To adequately (not a very precise term) a variogram from a sample variogram you need a moderate number of half squared differences in each distance class particularly for the shorter distances. This means a compromise between the number of plotted points (i.e. distance classes) and the number of pairs of data locations used to determine each plotted point. In addition the spatial pattern of the data locations greatly affects the number of pairs in each distance class. How widely spread are the clusters from each other and within the clusters how widely spread are the three data locations. With your spatial pattern you will get a moderate number of pairs for a short lag distance and a somewhat larger number of pairs for a long distance but very little in between hence you have little information to determine the shape of the variogram.

When you say the variance doesn't change with distance either you are mixing up the terminology/ideas or misunderstanding what a variogram is. If the sample variogram appears to level off and be asymptotic to a horizontal line then the height of the line is the estimated value of the variance of the implicit random function. To use a variogram the implicit random function must be intrinsic stationary, this is not a property of data. If you are referring to the variation of the half squared differences that is quite different.

If in fact the graph of the sample variogram looks like a horizontal line then as noted by others then this would correspond to a pure nugget effect and the best interpolator would simply be the arithmetic average of the data (no spatial correlation)

"Random trend" is nonsensical. If the data appears to have trend then you must use the residuals (from a trend surface fit) to compute the sample variogram. If there is a trend in the data then the sample variogram will appear to increase quadratically and there are no variogram models that increase that rapidly

The thin plate spline is in fact the same as ordinary kriging, it just looks different.

There are a few examples in the  geostatistics literature with data sets as small as this but in most cases the investigator has other information/knowledge about the phenomenon that generated the data

The local estimator method with the third-order polynomial function is the best method.