krigeST exploits the separability of S/T covariances to make kriging
more efficient, but does not provide simulation - it would be an option,
and is maybe possible for Gaussian simulation (as you're solving the
problem globally) for data sets of limited size.
For indicator kriging, afaics the algorithm needs to be sequential,
which is not implemented.
predict() in gstat allows for 3D model, meaning for S/T data you can use
x/y/t coordinates with the metric model; it uses the sequential
algorithm, so can do Gaussian as well as indicator simulation. It is
advisable to rescale the t dimension such that correlation range in x/y
is simular to that in t direction.