The ability of short pulse laser ablation in liquids to produce clean colloidal nanoparticles and unusual surface morphology has been employed in a broad range of practical applications. In this paper, we report the results of large-scale molecular dynamics simulations aimed at revealing the key processes that control the surface morphology and nanoparticle size distributions by pulsed laser ablation in liquids. The simulations of bulk Ag targets irradiated in water are performed with an advanced computational model combining a coarse-grained representation of liquid environment and an atomistic description of laser interaction with metal targets. For the irradiation conditions that correspond to the spallation regime in vacuum, the simulations predict that the water environment can prevent the complete separation of the spalled layer from the target, leading to the formation of large subsurface voids stabilized by rapid cooling and solidification. The subsequent irradiation of the laser-modified surface is found to result in a more efficient ablation and nanoparticle generation, thus suggesting the possibility of the incubation effect in multipulse laser ablation in liquids. The simulations performed at higher laser fluences that correspond to the phase explosion regime in vacuum reveal the accumulation of the ablation plume at the interface with the water environment and the formation of a hot metal layer. The water in contact with the metal layer is brought to the supercritical state and provides an environment suitable for nucleation and growth of small metal nanoparticles from metal atoms emitted from the hot metal layer. The metal layer itself has limited stability and can readily disintegrate into large (tens of nanometers) nanoparticles. The layer disintegration is facilitated by the Rayleigh-Taylor instability of the interface between the higher density metal layer decelerated by the pressure from the lighter supercritical water. The nanoparticles emerging from the layer disintegration are rapidly cooled and solidified due to the interaction with water environment, with a cooling rate of ∼2 × 1012 K/s observed in the simulations. The computational prediction of two distinct mechanisms of nanoparticle formation yielding nanoparticles with different characteristic sizes provides a plausible explanation for the experimental observations of bimodal nanoparticle size distributions in laser ablation in liquids. The ultrahigh cooling and solidification rates suggest the possibility for generation of nanoparticles featuring metastable phases and highly nonequilibrium structures.