I don't believe this is possible. Here's a sketch of a proof. Consider a level set $\Sigma_c$ of $f=c$. The gradient $\nabla f$ is a constant length vector perpendicular to $\Sigma_c$ at each point of $\Sigma_c$ (let's assume $|\nabla f|=1$). Also, the vector field $\nabla f$ is divergence-free, since $f$ is harmonic. Now, flow by the vector field $\nabla f$ for time $t$ takes $\Sigma_c$ into $\Sigma_{c+t}$. Since the vectors are constant length, this gives a local orthogonal coordinate system about $\Sigma_c$, which therefore is Fermi coordinates. So the flow lines of $\nabla f$ are geodesics by the Gauss lemma. This flow is also volume preserving since $\nabla f$ is harmonic, and therefore also preserves the area of $\Sigma_{c+t}$ since $|\nabla f|=1$. One then computes that the principal curvatures of $\Sigma_c$ at each point are $\pm 1$. Thus, by Gauss' equation, $\Sigma_c$ is isometric to $\mathbb{H}^2_{-2}$. However, by a result of Doug Moore (generalizing a classic result of Hilbert), there is no isometric immersion of $\mathbb{H}^2_{-2}$ into $\mathbb{H}^3_{-1}$.
It's possible that I've made an error here though, since I haven't checked the computations.