Spatio-temporal block kriging with R package gstat?

In sp, SpatialPoints*, SpatialPixels* and SpatialGrid* (with * omitted or replaced by DataFrame) do support more than 2 spatial dimensions, as OP has done, but SpatialPolygons* and SpatialLines* do not. With gstat you can do 3-D block kriging with 3-D blocks (using block = c(10,10,10)), but you cannot do this for non-rectangular blocks, as OP wants. It is perfectly OK to substitute time for the third dimension, but you are constrained to the metric ST variogram.

library(gstat)
vignette("st")

gives you more options for variogram models, but not for predicting block mean values (this is FYI, not an answer to the question).

The only answer to the question would be to do 3D conditional simulations, and aggregate point values over your arbitrary 3D (2D polygon + time extent) blocks. Tedious, but possible; also only along the 3D path, not along the path described in the ST vignette (krigeST does not do simulation - yet!).

You really need to do a bit more research on your methodology and read the documentation to understand the structure of sp S4 class objects and interaction of sp objects with relevant gstat functions. In the sp Vignette there is a detailed explanation of the difference between SpatialPolygons (only polygon topology) and SpatialPolygonDataFrame (polygons with attributes) objects.

What you are explaining is not block Kriging and using time as an attribute does not result in a spatial-temporal estimate. The conceptual methodology you describe is quite invalid. Using polygons or polygon centroids violates the Kriging assumptions of a uniform random field, anisotropy and nonstationarity.

Here is a nice gstat vignette on spatial-temporal models using the interface to the spacetime package. I should also note that the constrainedKriging package can conduct block Kriging on arbitrary shaped blocks using a nonstationary mean function and an isotropic weakly stationary variogram.

That said, to answer your question, you can pass an sp SpatialPointsDataFrame object directly to a variogram/Kriging model in gstat. In this type of sp object, attributes reside in the "data" slot and are already attached to the coordinates via the internal S4 class structure.

# COERCE meuse DATAFRAME TO sp SpatialPointsDataFrame OBJECT
require(gstat)
data(meuse)
coordinates(meuse) <- ~ x + y
head(meuse@data)

# CREATE SEMIVARIOGRAM USING THE zinc ATTRIBUTE
# NOTE: THERE IS NO ARGUMENT FOR A "4th DIM"
v <- variogram(log(zinc) ~ 1, meuse)
plot(v, type = "l")