A search is performed for a scalar resonance X decaying into a Standard Model Higgs boson H and an additional scalar S, each decaying into a bottom-antibottom-quark pair, giving a $b\bar{b}b\bar{b}$ final state. Two event topologies are considered, depending on whether the Higgs boson is reconstructed from a single large-radius jet or from a pair of small-radius jets. The analysis uses proton-proton collision data at $\sqrt{s}=13$ TeV, collected by the ATLAS experiment at the Large Hadron Collider during Run 2, corresponding to an integrated luminosity of 126-136 fb${}^{-1}$ depending on the topology and the associated triggers. The search covers masses of 300 GeV-6 TeV for X and 70 GeV-5 TeV for S. The dominant multijet and $t\bar{t}$ backgrounds are estimated with a hierarchical model using Gaussian processes and normalizing flows. No significant excess over the Standard Model background expectation is observed. Upper limits at 95% confidence level are set on the production cross-section times branching ratio, ranging from 0.7 fb to 2.6 pb.