We developed a numerical model of compression of asperities by pressure solution (PS). The dissolution rate along the contact was determined by (1) computing the normal stress distribution from the present shape of the asperities, and (2) solving the diffusion equation inside the fluid-saturated solid-solid interface, including local dissolution source terms corresponding to the stress field previously determined. The change in shape of the asperities during an infinitesimal time interval can then be calculated and the entire procedure repeated as many times as desired. We simulated PS compaction of axisymmetric asperities with different sizes and shapes under various temperatures and loads, and using different values of the interface diffusion coefficient. Our results show that as the contact flattens and grows during PS, the initial elastic deformation is partially relaxed and the stress transferred from the contact centre to the edge. Transient stress release remains significant during an extended period and could strongly influence the interpretation of laboratory experiments. In our model, dissolution and interface diffusion are not sequentially combined as is usually assumed, making it impossible to identify a single rate-limiting process. The convergence rate of the asperities was approximately proportional to the mean effective stress at the contact, and depended in a complex way on the contact radius, but was not sensitive to the asperity size. We also estimated the conditions under which undercutting at the contact edge is triggered.