Kriging with gstat: including a random effect

I am getting acquainted with gstat and kriging more generally in R. The end goal is to interpolate temperature across my spatial grid, and I want to be able to do so for different time points (average daily temperature over several years). So, I've been imagining an interpolation that estimates temperature and considers something like the season or year as a random effect. I've been exploring universal kriging, but I'm struggling to understand the formula set up.

I've included some example code using the meuse data (from sp). When I assign a variogram, I can state that my dependent variable is predicted by other variables and include random effects (rather than using something like zinc~1).

However, when I do this with krige(), the only alternative formula to zinc~1 seems to be zinc~x+y. I believe the formula zinc~x+y is what allows anisotropy and avoids the stationarity assumption, but why can't I include other data here like I can with the variogram? Is it because it is already accounted for in the variogram, or am I not understanding something about kriging more generally?

library(sp)
library(gstat)

data(meuse)
coordinates(meuse) <- ~ x+y
data(meuse.grid)
coordinates(meuse.grid) <- ~ x+y

# assigning and fitting the variogram
sample.vgm <- variogram(zinc~copper+(1|dist.m), # the formula - copper is a fixed effect, dist.m is a random effect 
                     meuse) # SPDF dataset

vgm.fit <- fit.variogram(sample.vgm, # the sample variogram 
                         model=vgm(model="Sph")) # assigning a model type

# kriging
z.krige <- krige(zinc~x+y, # the  formula
                 meuse, # the SPDF that has data
                 meuse.grid, # the SPDF to krige over
                 model=vgm.fit) 

# this doesn't work if the formula is anything besides ~1 or ~x+y
z.krige2 <- krige(zinc~copper+(1|dist.m),
                 meuse, 
                 meuse.grid, 
                 model=vgm.fit)

I am using R version 4.1.3 on Windows 10 (x64).

Solution

Based on this post, the explanatory variables apparently need to be loaded into the grid for Kriging to find them.

dist should work since it is already in the meuse grid:

# dist is in the grid
z.krige2 <- krige(zinc~dist,
  meuse, 
  meuse.grid, 
  model=vgm.fit)

For copper and dist.m, I've just added random values here:

# add random values for copper and dist.m into the grid
meuse.grid$copper <- runif(nrow(meuse.grid))
meuse.grid$dist.m <- runif(nrow(meuse.grid))

# now copper and dist.m variables are available
z.krige2 <- krige(zinc~copper+(1|dist.m),
  meuse, 
  meuse.grid, 
  model=vgm.fit)