Merging neutron stars are expected to produce hot, metastable remnants in rapid differential rotation, which subsequently cool and evolve into rigidly rotating neutron stars or collapse to black holes. Studying this metastable phase and its further evolution is essential for the prediction and interpretation of the electromagnetic, neutrino, and gravitational signals from such a merger. In this work, we model binary neutron star merger remnants and propose new rotation and thermal laws that describe postmerger remnants. Our framework is capable to reproduce quasiequilibrium configurations for generic equations of state, rotation and temperature profiles, including nonbarotropic ones. We demonstrate that our results are in agreement with numerical relativity simulations concerning bulk remnant properties like the mass, angular momentum, and the formation of a massive accretion disk. Because of the low computational cost for our axisymmetric code compared to full 3 + 1-dimensional simulations, we can perform an extensive exploration of the binary neutron star remnant parameter space studying several hundred thousand configurations for different equations of state.