This paper presents a Lagrangian stochastic approach for modeling field scale contaminant transport in the heterogeneous unsaturated geological environments. We derive the probability distribution function of the travel times of tracers and reactive solutes, without limiting it to be either normal or log-normal, and show how it can be used for computing the moments of solute fluxes and breakthrough curves. We also provide closed-form expressions for the travel time moments of tracers and reactive solutes characterized by linear equilibrium and non-equilibrium sorption kinetics. These expressions are useful for prediction and for interpretation of field experiments. In our derivation we account for soil heterogeneity by modeling the soil parameters, such as the saturated hydraulic conductivity and the pore size distribution parameter, as weakly stationary random space functions. Unlike previous studies, we account for the effect of the spatial variability of the water content and refrain from assuming that the saturation is practically constant in the unit gradient flow zone. The derivation is done using a first-order perturbative expansion of the flow equation, and is thus limited to small variability of the input parameters. The effects of spatial variability on the displacement and travel time statistics are discussed and demonstrated.