Simulations have been performed with a fairly narrow band numerical gravity wave model (higher-order nls type) and a computational domain of dimensions $128\times 128$ typical wavelengths. The simulations are initiated with $\thicksim6\times10^{4}$ fourier modes corresponding to truncated jonswap spectra and different angular distributions giving both short- and long-crested waves. A development of the spectra on the so-called benjamin–feir timescale is seen, similar to the one reported by dysthe et al. (J. Fluid Mech. vol. 478, 2003, p.1). The probability distributions of surface elevation and crest height are found to fit theoretical distributions found by tayfun (J. Geophys. Res. Vol. 85, 1980, p. 1548) very well for elevations up to four standard deviations (for realistic angular spectral distributions). moreover, in this range of the distributions, the influence of the spectral evolution seems insignificant. for the extreme parts of the distributions a significant correlation with the spectral change can be seen for very long-crested waves. For this case we find that the density of large waves increases during spectral change, in agreement with a recent experimental study by onorato et al. (J. Fluid Mech. 2004 submitted).