Abstract.We investigate the critical behavior of the three-dimensional random-field Ising model (RFIM) with a Gaussian field distribution at zero temperature. By implementing a computational approach that maps the ground-state of the RFIM to the maximum-flow optimization problem of a network, we simulate large ensembles of disorder realizations of the model for a broad range of values of the disorder strength h and system sizes V = L 3 , with L ≤ 156. Our averaging procedure outcomes previous studies of the model, increasing the sampling of ground states by a factor of 10 3 . Using well-established finitesize scaling schemes, the fourth-order's Binder cumulant, and the sample-to-sample fluctuations of various thermodynamic quantities, we provide high-accuracy estimates for the critical field hc, as well as the critical exponents ν, β/ν, andγ/ν of the correlation length, order parameter, and disconnected susceptibility, respectively. Moreover, using properly defined noise to signal ratios, we depict the variation of the self-averaging property of the model, by crossing the phase boundary into the ordered phase. Finally, we discuss the controversial issue of the specific heat based on a scaling analysis of the bond energy, providing evidence that its critical exponent α ≈ 0 − .