First-order analysis, based on a stochastic continuum presentation of flow and a general Lagrangian description of transport, was used to investigate the effects of several characteristics of a bimodal, spatially heterogeneous, variably saturated formation, on solute breakthrough curves, under steady-state, gravity-dominated unsaturated flow conditions. The bimodal formation was viewed as a mixture of two populations (background soil and embedded soil) of differing hydraulic properties and differing spatial structures. Results of the present analysis, restricted to relatively small inclusions' volume fraction, suggest that as in saturated flow, for a formation of given statistics and mean water saturation, both the spreading of the mean solute breakthrough curve and the uncertainty in its prediction may be appreciable, especially in bimodal formations in which: (i) the contrast between the mean properties of the two subdomains is relatively large, (ii) the inclusions' volume fraction is relatively large, and (iii) the characteristic length scales of the inclusions are relatively large as compared with those of the heterogeneity of the background soil. Results of this study, unique to unsaturated flow conditions, stem from the inherent concave nature of the log-conductivity variance–mean pressure head relationships in bimodal, variably saturated formations. Consequently, under unsaturated flow conditions, both the spreading of the mean solute discharge and the uncertainty in its prediction are concave functions of mean saturation. Starting with saturated formation, they initially decrease with decreasing mean saturation, reach their minima, and then increase as mean saturation further decreases, exceeding their counterparts in saturated flow conditions. Implications of the results with respect to the problem of groundwater contamination where risk assessment is required, are briefly discussed.