Let \(b^p\), \(1\leq p< \infty\), be the harmonic Bergman space in the unit ball of \(\mathbb{R}^n\). The following bounds for (derivatives of) the harmonic Bergman kernel \(R(x,y)\) are obtained: \[ |\partial_x^\alpha \partial_y^\alpha R(x,y)|\leq C(1-2x\cdot y+|x|^2|y|^2)^{-(n+|\alpha|+ |\beta|)/2}. \] They are applied, together with some reproducing formulas, to solve Gleason's problem in \(b^p\) and in the harmonic Bloch space, and to get a characterization of these spaces in terms of derivative norms.