Given a linear partial differential equation (PDE), I will discuss the design of Gaussian process (GP) priors for the approximation of its solutions as well as its physical parameters. As standard PDE theory is based on weak or distributional formulations, I will describe necessary and sufficient conditions on the kernel of a centered GP so that its realisations solve the PDE almost surely, in the distributional sense. I will then describe necessary and sufficient conditions so that the samples of a centered GP lie in a given Sobolev space, as the latter are particularly well-suited for the study of PDEs. Importantly, both results do not make any continuity assumptions on the GP model, as PDE theory shows that such assumptions are sometimes too restrictive. If time permits, I will describe an application of such GP models for the approximation of the initial data of the 3D wave equation, as is e.g. the goal of photo-acoustic tomography.