We initiate a detailed study of two-parameter Besov spaces on the unit ball of ℝn consisting of harmonic functions whose sufficiently high-order radial derivatives lie in harmonic Bergman spaces. We compute the reproducing kernels of those Besov spaces that are Hilbert spaces. The kernels are weighted infinite sums of zonal harmonics and natural radial fractional derivatives of the Poisson kernel. Estimates of the growth of kernels lead to characterization of integral transformations on Lebesgue classes. The transformations allow us to conclude that the order of the radial derivative is not a characteristic of a Besov space as long as it is above a certain threshold. Using kernels, we define generalized Bergman projections and characterize those that are bounded from Lebesgue classes onto Besov spaces. The projections provide integral representations for the functions in these spaces and also lead to characterizations of the functions in the spaces using partial derivatives. Several other applications follow from the integral representations such as atomic decomposition, growth at the boundary and of Fourier coefficients, inclusions among them, duality and interpolation relations, and a solution to the Gleason problem. © 2016 World Scientific Publishing Company.