Thanks, I like the idea of treating the vector field as the gradient of a 3D function. Just to clarify on step 3: if I understand correctly what you describe is an optimization over the weights assigned to the basis vectors, and the output would be a single field that most closely matches my observations. So to find the conditional distribution would I need to use something like Laplace's method, or did you have something else in mind?

Thanks for your comment, I've added some further details. Either way, I'd be interested to know which interpretation would make the problem straightforward.