  A search for pair production of the supersymmetric partners of the Higgs
boson (higgsinos $\tilde{H}$) in gauge-mediated scenarios is reported. Each
higgsino is assumed to decay to a Higgs boson and a gravitino. Two
complementary analyses, targeting high- and low-mass signals, are performed to
maximize sensitivity. The two analyses utilize LHC $pp$ collision data at a
center-of-mass energy $\sqrt{s} = 13$ TeV, the former with an integrated
luminosity of 36.1 fb$^{-1}$ and the latter with 24.3 fb$^{-1}$, collected with
the ATLAS detector in 2015 and 2016. The search is performed in events
containing missing transverse momentum and several energetic jets, at least
three of which must be identified as $b$-quark jets. No significant excess is
found above the predicted background. Limits on the cross-section are set as a
function of the mass of the $\tilde{H}$ in simplified models assuming
production via mass-degenerate higgsinos decaying to a Higgs boson and a
gravitino. Higgsinos with masses between 130 and 230 GeV and between 290 and
880 GeV are excluded at the 95% confidence level. Interpretations of the limits
in terms of the branching ratio of the higgsino to a $Z$ boson or a Higgs boson
are also presented, and a 45% branching ratio to a Higgs boson is excluded for
$m_{\tilde{H}} \approx 400$ GeV.
