Further to my answer in the comments, if you wanted to look at the variance of the ice curves, you could bootstrap them like this:
library(pdp)
library(randomForest)
library(ICEbox)
data(boston)
X <- as.data.frame(model.matrix(cmedv ~ ., data=boston)[,-1])
y <- model.response(model.frame(cmedv ~ ., data=boston))
boston.rf <- randomForest(x=X, y=y)
bice <- ice(boston.rf, X=X, predictor = "lstat") 
res <- NULL
for(i in 1:1000){
  inds <- sample(1:nrow(bice$ice_curves), 
                 nrow(bice$ice_curves), 
                 replace=TRUE)
  res <- rbind(res, colMeans(bice$ice_curve[inds, ]))
}
out <- data.frame(
  fit = colMeans(bice$ice_curves), 
  lwr = apply(res, 2, quantile, .025),
  upr = apply(res, 2, quantile, .975), 
  x=bice$gridpts
)
library(ggplot2)
ggplot(out, aes(x=x, y=fit, ymin=lwr, ymax=upr)) + 
  geom_ribbon(alpha=.25) + 
  geom_line() + 
  theme_bw() + 
  labs(x="lstat", y="Prediction")

Or, you could look at the different quantiles of the ice plots for each evaluation point.
tmp <- t(apply(bice$ice_curves, 
             2, 
             quantile, c(0, .025, .05, .1, .25, .5, .75, .9, .95, .975, 1)))
head(tmp)
tmp <- as.data.frame(tmp)
names(tmp) <- c("l1", "l2", "l3", "l4", "l5", 
                "med", "u1", "u2", "u3", "u4", "u5")
tmp$x <- bice$gridpts
ggplot(tmp, aes(x=x, y=med)) + 
  geom_ribbon(aes(ymin=l1, ymax=u1), alpha=.2) + 
  geom_ribbon(aes(ymin=l2, ymax=u2), alpha=.2) + 
  geom_ribbon(aes(ymin=l3, ymax=u3), alpha=.2) + 
  geom_ribbon(aes(ymin=l4, ymax=u4), alpha=.2) + 
  geom_ribbon(aes(ymin=l5, ymax=u5), alpha=.2) + 
  geom_line() + 
  theme_bw() + 
  labs(x="lstat", y="Prediction")
