I actually just stumbled on this paper and got very interested in this. Since this was posted 2 years ago, I am not sure OP is still looking for an explanation of how this works.
The idea is that it is essentially a multi-output Gaussian process, except that the full covariance function is specified once you specify the covariance function for the function observations. This is because the derivative is a linear operation and the resulting process for derivative is also Gaussian. Furthermore, the relationship between function and derivative imposes constraints on what the covariance function looks like. So if you have a function with input $\mathbf{x}=(x_1, x_2,...,x_p)$ (each one of your input is $p$-dimensional) , you are then fitting a Gaussian process with $(p+1)$ outputs ($1$ for function values and $p$ partial derivatives) and the covariance function is as specified in the paper.
The problem with adding derivative information is drastic increase in computational time if your function has high dimensional input. When you only have function observation, you only need to solve for a $n\times n$ matrix assuming you have $n$ data points. But with $p$-dimensional input, you now need to solve for a $n(p+1)\times n(p+1)$ matrix which will be extremely costly for large $p$. So it depends on whether you're willing to have this tradeoff between computational cost and accuracy.
Hope this helps.