  We model the local stellar velocity field using position and velocity
measurements for 4M stars from the second data release of Gaia. We determine
the components of the mean or bulk velocity in ~27k spatially-defined bins. Our
assumption is that these quantities constitute a Gaussian process where the
correlation between the bulk velocity at different locations is described by a
simple covariance function or kernel. We use a sparse Gaussian process
algorithm based on inducing points to construct a non-parametric, smooth, and
differentiable model for the underlying velocity field. We estimate the Oort
constants A, B, C, and K and find values in excellent agreement with previous
results. Maps of the velocity field within 2 kpc of the Sun reveal complicated
substructures, which provide clear evidence that the local disk is in a state
of disequilibrium. We present the first 3D map of the divergence of the stellar
velocity field and identify regions of the disk that may be undergoing
compression and rarefaction.
