The hadronic production of a Higgs boson (H) in association with b jets will play an important role in investigating the Higgs-boson couplings to Standard Model particles during Run II of the CERN Large Hadron Collider, and could in particular reveal the presence of anomalies in the assumed hierarchy of Yukawa couplings to the third-generation quarks. A very high degree of accuracy in the theoretical description of this process is crucial to implement the rich physics program that could lead to either direct or indirect evidence of new physics from Higgs-boson measurements. Aiming for accuracy in the theoretical modeling of H+b-jet production, we have interfaced the analytic Next-to-Leading-Order QCD calculation of H-bottom-antibottom production with parton-shower Monte Carlo event generators in the POWHEG BOX framework. In this paper we describe the most relevant aspects of the implementation and present results for the production of H+1 b jet, H+2 b jets, and $H$ with no tagged b jets, in the form of kinematic distributions of the Higgs boson, of the b jets, and of the non-b jets, at the 13 TeV Large Hadron Collider. The corresponding code is part of the public release of the POWHEG BOX.