This paper presents a search for pair production of higgsinos, the supersymmetric partners of the Higgs bosons, in scenarios with gauge-mediated supersymmetry breaking. Each higgsino is assumed to decay into a Higgs boson and a nearly massless gravitino. The search targets events where each Higgs boson decays into$b\overline{b}$, leading to a reconstructed final state with at least three energetic$b$-jets and missing transverse momentum. Two complementary analysis channels are used, with each channel specifically targeting either low or high values of the higgsino mass. The low-mass (high-mass) channel exploits$126(139){\mathrm{fb}}^{-1}$of$\sqrt{s}=13\mathrm{TeV}$data collected by the ATLAS detector during Run 2 of the Large Hadron Collider. No significant excess above the Standard Model prediction is found. At 95% confidence level, masses between 130 GeV and 940 GeV are excluded for higgsinos decaying exclusively into Higgs bosons and gravitinos. Exclusion limits as a function of the higgsino decay branching ratio to a Higgs boson are also reported.