We report on the investigation of height distributions (HDs) and spatial covariances of two-dimensional surfaces obtained from extensive numerical simulations of the celebrated Clarke-Vvedensky (CV) model for homoepitaxial thin film growth. In this model, the effect of temperature, deposition flux, and strengths of atom-atom interactions are encoded in two parameters: the diffusion to deposition ratio R = D/F and ε, which is related to the probability of an adatom "breaking" a lateral bond. We demonstrate that the HDs present a strong dependence on both R and ε, and even after the deposition of 10 5 monolayers (MLs) they are still far from the asymptotics in some cases. For instance, the temporal evolution of the HDs' skewness (kurtosis) displays a pronounced minimum (maximum), for small R and ε, and only at long times it passes to increase (decrease) toward its asymptotic value. However, it is hard to determine whether they converge to a single value or different nonuniversal ones. For large R and/or ε, on the other hand, these quantities clearly converge to the values expected for the Villain-Lai-Das Sarma (VLDS) universality class. A similar behavior is observed in the spatial covariances, but with weaker finite-time effects, so that rescaled curves of them collapse quite well with the one for the VLDS class at long times. Simulations of a model with limited mobility of particles, which captures some essential features of the CV model in the limit of irreversible aggregation (ε = 0), reveal a similar scenario. Overall, these results point out that the study of fluctuations in homoepitaxial thin films' surfaces can be a very difficult task and shall be performed very carefully, once typical experimental films have 10 4 MLs, so that their HDs and covariances can be in the realm of transient regimes.