The aim of this paper is to analyze the statistical properties of solute concentration in natural aquifers as sampled in observation wells, having a small diameter in comparison with the characteristic size of the heterogeneity in hydraulic properties. The analysis, in Langragian framework, takes advantage of the "reverse formulation", where, instead of considering the destination of the injected particles, the origin of the particle being sampled is sought. In the case of small values of the log-conductivity variance sigma(Y)(2) it allows the derivation of an analytical expression for concentration mean, variance and pdf, while for aquifer characterized by high value in sigma(Y)(2), a numerical analysis based on a Monte Carlo approach using a reverse scheme is developed and applied for values of sigma(Y)(2) up to 2. In this case, the use of a Beta function to fit the concentration pdf proves valid for practical applications. The comparison between the numerical and the analytical results defines the range of validity of the analytical ones. The relative role of large-scale dispersion processes and pore-scale effects is analyzed in terms of global variance in order to point out limits and accuracy of the Eulerian scheme in comparison with the Lagrangian one.