It is common to represent solute tranport in heterogeneous formations in terms of the resident concentration C (x, t), regarded as a random space function. The present study investigates the alternative representation by q , the solute mass flux at a point of a control plane normal to the mean flow. This representation is appropriate for many field applications in which the variable of interest is the mass of solute discharged through a control surface. A general framework to compute the statistical moments of q and of the associated total solute discharge Q and mass M is established. With x the direction of the mean flow, a solute particle is crossing the control plane at y = η, z = ζ and at the travel (arrival) time τ. The associated expected solute flux value is proportional to the joint probability density function (pdf) g ₁ of η, ζ and τ, whereas the variance of q is shown to depend on the joint pdf g ₂ of the same variables for two particles. In turn, the statistical moments of η, ζ and τ depend on those of the velocity components through a system of stochastic ordinary differential equations. For a steady velocity field and neglecting the effect of pore‐scale dispersion, a major simplification of the problem results in the independence of the random variables η, ζ and τ. As a consequence, the pdf of η and ζ can be derived independently of τ. A few approximate approaches to derive the statistical moments of η, ζ and τ are outlined. These methods will be explored in paper 2 in order to effectively derive the variances of the total solute discharge and mass, while paper 3 will deal with the nonlinear effect of the velocity variance upon the moments of η, ζ and τ