OK, I think I have a way to make an approximation that holds up pretty well by using Laplace's method twice.
Begin by using the representation
$$K_0(u) = \int_0^{\infty} dt \, e^{-u\cosh{t}}$$
Then the integral you seek may be written, when we reverse the order of integration, as
$$\int_0^{\infty} dk \, K_0\left(\sqrt{a (k^2+b)}\right) = \int_0^{\infty} dt \, \int_0^{\infty} dk \, e^{-\sqrt{a} \cosh{t} \sqrt{k^2+b}}$$
This is justified in that both integrals clearly converge absolutely. Now, to apply Laplace's method on the integral over $k$, we may assume that $a \cosh{t}$ is sufficiently large, and we use the approximation that
$$\sqrt{k^2+b} \sim \sqrt{b} + \frac{k^2}{2 \sqrt{b}}$$
We end up with a familiar integral that may be evaluated immediately, and we are left with a single integral:
$$\int_0^{\infty} dk \, K_0\left(\sqrt{a (k^2+b)}\right) \sim \sqrt{\frac{\pi \sqrt{b}}{2 \sqrt{a}}} \int_0^{\infty} \frac{dt}{\sqrt{\cosh{t}}} e^{-\sqrt{a b} \cosh{t}} \quad (a \to \infty)$$
At this point I should say that I am not providing an error estimate, although deriving one should be fairly straightforward. (This involves providing an additional term of the Taylor expansion of the square root above, and then Taylor expanding the exponential in the integrand.)
Now we are left with the problem of evaluating the above integral, which looks only a little less problematic than the original one. Nevertheless, we may apply Laplace's method again, as we are claiming that $a$ is sufficiently large, and $b \gt 0$. In this case, we use the Taylor expansion of $\cosh{t} \sim 1+t^2/2$ and we again end up with a familiar integral. The final result is
$$\int_0^{\infty} dk \, K_0\left(\sqrt{a (k^2+b)}\right) \sim \frac{\pi}{2 \sqrt{a}} e^{-\sqrt{a b}} \quad (a \to \infty)$$
Some random numerical samples - even with moderate values of $a=1$ and $b=1$ in WA have verified the correctness of this approximation.