How do I take groundwater flow direction into consideration when 3D kriging?

Kriging can’t specifically take flow direction and magnitude as terms, but you can account for flow velocity (velocity is a vector and therefore includes speed and direction) using a fully anisotropic variogram.

The Advanced Variography in krig_3d class specifically addresses this complex topic.

One of the limitations is that kriging can’t have a preferential “downstream” direction.  This is rarely ever an issue, but you need to be aware that a downstream spill would interpolate up-flow just as preferentially as down-flow since the manner by which kriging assigns preferential directions uses symmetrical ellipsoids.