Despite its ubiquitous character and relevance in many branches of science and engineering, nucleation from solution remains elusive. In this framework, molecular simulations represent a powerful tool to provide insight into nucleation at the molecular scale. In this work, we combine theory and molecular simulations to describe urea nucleation from aqueous solution. Taking advantage of well-tempered metadynamics, we compute the freeenergy change associated to the phase transition. We find that such a free-energy profile is characterized by significant finite-size effects that can, however, be accounted for. The description of the nucleation process emerging from our analysis differs from classical nucleation theory. Nucleation of crystal-like clusters is in fact preceded by large concentration fluctuations, indicating a predominant two-step process, whereby embryonic crystal nuclei emerge from dense, disordered urea clusters. Furthermore, in the early stages of nucleation, two different polymorphs are seen to compete.T he nucleation of crystals from solution plays an important role in a variety of chemical and engineering processes. However, the very small scale that characterizes the early stages of nucleation makes its experimental study rather challenging. In this regard, molecular simulations play an important role and much work has been devoted to the study of homogeneous nucleation in simple model systems like Lennard-Jones particles or hard spheres (1-3). More recently, growing attention has been paid to the computer simulation of nucleation from solution of organic and inorganic materials (4-11).However, these simulations have to face an important challenge. Nucleation is a typical example of a rare event occurring on a timescale that is much longer than what atomistic simulations can typically afford. The problem is most easily understood if we focus on molecular dynamics (MD), which is the sampling method used here but it plagues other methods as well. In MD, the shortest timescale to be used for a correct integration of the equation of motion is dictated by the fastest atomic motions such as vibrations or rotations. Due to this constraint in the integration step, present MD simulations can reach timescales up to microseconds unless specific hardware, accessible only to few, is used. Even in this case the scale of the accessible times can be stretched only to the milliseconds range, a far cry from the much longer scale of nucleation phenomena.These timescale limitations affect molecular simulations in several other research fields, such as ligand protein binding, protein folding, or slow chemical reactions, to name just a few. To overcome this limit, many enhanced sampling methods have been proposed (12-21). Several of them are based on the application of a suitable external bias potential (12-14) that speeds up configurational sampling and permits free energies to be evaluated and transition rates to be computed (22)(23)(24). Here, we shall use well-tempered (WT) metadynamics to enhance the nucleation...