  We present up-to-date matched predictions for the $b\bar{b}H$ inclusive cross
section at the LHC at $\sqrt{s}=13$ TeV. Using a previously developed method,
our predictions consistently combine the complete NLO contributions that are
present in the 4-flavor scheme calculation, including finite b-quark mass
effects as well as top-loop induced $Y_b Y_t$ interference contributions, with
the resummation of collinear logarithms of $m_b/m_H$ as present in the 5-flavor
scheme calculation up to NNLO. We provide a detailed estimate of the
perturbative uncertainties of the matched result by examining its dependence on
the factorization and renormalization scales, the scale of the Yukawa coupling,
and also the low b-quark matching scale in the PDFs. We motivate the use of a
central renormalization scale of $m_H$/2, which is halfway between the values
typically chosen in the 4-flavor and 5-flavor scheme calculations. We evaluate
the parametric uncertainties due to the PDFs and the b-quark mass, and in
particular discuss how to systematically disentangle the parametric $m_b$
dependence and the unphysical b-quark matching scale dependence. Our best
prediction for the $b\bar{b} H$ production cross section in the Standard Model
at 13 TeV and for $m_H$ = 125 GeV is $\sigma(b\bar b H) = 0.52 \,{\rm pb}\, [1
\pm 9.6\% {\rm (perturbative)}\,{}^{+2.9\%}_{-3.6\%} {\rm (parametric)}]$. We
also provide predictions for a range of Higgs masses $m_H\in [50, 750]$ GeV.
Our method to compute the matched prediction and to evaluate its uncertainty
can be readily applied to other heavy-quark-initiated processes at the LHC.
