A search for Higgs bosons that decay into a b quark-antiquark pair and are accompanied by at least one additional b quark is performed with the CMS detector. The data analyzed were recorded in proton-proton collisions at a centre-of-mass energy of 13~\mathrm{TeV} at the LHC, corresponding to an integrated luminosity of 35.7~\mathrm{fb}^{-1}. The final state considered is particularly sensitive to signatures of a Higgs sector beyond the standard model, as predicted by the minimal supersymmetric standard model (MSSM) and the two Higgs doublet model (2HDM) with large values of the parameter . No signal above the standard model background expectation is observed. Stringent upper limits on the cross section times branching fraction are set for Higgs bosons with masses up to 1300~\mathrm{GeV} at confidence level. The results are interpreted within several MSSM and 2HDM scenarios. In the hMSSM scenario, upper limits on are obtained, ranging from to for Higgs masses from to 900~\mathrm{GeV}. In the flipped 2HDM scenario, similar upper limits on are set over the full range and for Higgs masses from to 850~\mathrm{GeV}.