© 2018 CERN for the ATLAS Collaboration. Published by the American Physical Society under the terms of the »https://creativecommons.org/licenses/by/4.0/» Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP 3 . A search for pair production of the supersymmetric partners of the Higgs boson (higgsinos 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 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 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 mH≈400 GeV.