In recent years, statistical theory has been used to compute the ensemble mean and variance of solute concentration in aquifer formations with second‐order stationary velocity fields. The merit of accurately estimating the mean and variance of concentration, however, remains unclear without knowing the shape of the probability density function (pdf). In a setup where a conservative solute is continuously injected into a domain, the concentration is bounded between zero and the concentration value in the injected solution. At small travel distances close to the fringe of the plume, an observation point may fall into the plume or outside, so that the statistical concentration distribution clusters at the two limiting values. Obviously, this results in non‐Gaussian pdf's of concentration. With increasing travel distance, the lateral plume boundaries are smoothed, resulting in increased probability of intermediate concentrations. Likewise, averaging the concentration in a larger sampling volume, as typically done in field measurements, leads to higher probabilities of intermediate concentrations. We present semianalytical results of concentration pdf's for measurements with point‐like or larger support volumes based on stochastic theory applied to stationary media. To this end, we employ a reversed auxiliary transport problem, in which we use analytical expressions for first and second central spatial lateral moments with an assumed Gaussian pdf for the uncertainty of the first lateral moment and Gauss‐like shapes in individual cross sections. The resulting concentration pdf can be reasonably fitted by beta distributions. The results are compared to Monte Carlo simulations of flow and steady state transport in 3‐D heterogeneous domains. In both methods the shape of the concentration pdf changes with distance to the contaminant source: Near the source, the distribution is multimodal, whereas it becomes a unimodal beta distribution far away from the contaminant source. The semianalytical and empirical pdf's differ slightly, which we contribute to the numerical artifacts in the Monte Carlo simulations but also to hard assumptions made in the semianalytical approach. Our results imply that geostatistical techniques for interpolation and other statistical inferences based on Gaussian distributions, such as kriging and cokriging, may be feasible only far away from the contaminant source. For calculations near the source, the beta‐like distribution of concentration should be accounted for.