Approximate analytic solutions are developed to predict the transport of inert solutes in variably saturated, randomly heterogeneous porous media subject to a random boundary flux condition. The mean and variance of solute mass flux past a series of vadose zone control planes are determined using a Lagrangian approach. The method requires the specification of an advective travel time probability density function (pdf) between the soil surface and the control planes and an impulse response function which quantifies the specific flow and transport processes that affect the solute along an individual streamline in the soil. In this study an advective-dispersive impulse response function for inert solute is assumed, and an advective travel time pdf is derived from the statistics of the boundary flux and soil properties using a one-dimensional transient vadose zone flow model <Foussereau et al., this issue>. Monte Carlo simulations are conducted using random boundary flux sequences and spatially correlated random hydraulic conductivity fields in a multidimensional variably saturated flow and transport code to test the analytic model assumptions. The analytic model accurately predicts the ensemble mean solute flux breakthrough curves for a variety of soil-climate combinations. The general behavior of the coefficient of variation of solute flux is also accurately reproduced; however, the minimum coefficient of variation at the mean center of mass was slightly underpredicted in most cases. Results show that for shallow soils in humid climates, boundary flux variability dominates soil variability in determining the uncertainty of solute transport predictions in the vadose zone. ¿ 2000 American Geophysical Union