We use a one-dimensional adaptive-grid thermal evolution code to model Kuiper belt objects Charon, Orcus and Salacia and compare their measured bulk densities with those resulting from evolutionary calculations at the end of 4.6 Gyr. Our model assumes an initial homogeneous composition of mixed ice and rock, and follows the multiphase flow of water through the porous rocky medium, consequent differentiation and aqueous chemical alterations in the rock. Heating sources include long-lived radionuclides, serpentinization reactions, release of gravitational potential energy due to compaction, and crystallization of amorphous ice. The density profile is calculated by assuming hydrostatic equilibrium to be maintained through changes in composition, pressure and temperature. To this purpose, we construct an equation of state suitable for porous icy bodies with radii of a few hundred km, based on the best available empirical studies of ice and rock compaction, and on comparisons with rock porosities in Earth analog and Solar System silicates. We show that the observed bulk densities can be reproduced by assuming the same set of initial and physical parameters, including the same rock/ice mass ratio for all three bodies. We conclude that the mass of the object uniquely determines the evolution of porosity, and thus explains the observed differences in bulk density. The final structure of all three objects is differentiated, with an inner rocky core, and outer ice-enriched mantle. The degree of differentiation, too, is determined by the object’s mass.