It is commonly assumed in stochastic solute (advective) transport models that either the velocity field is stationary (statistically homogeneous) or the mean flow is unidirectional. In this study, using a Lagrangian approach, we develop a general stochastic model for transport in variably saturated flow in randomly heterogeneous porous media. The mean flow in the model is multidirectional, and the velocity field can be nonstationary (with location-dependent statistics). The nonstationarity of the velocity field may be caused by statistical nonhomogeneity of medium properties or complex boundary configurations. The particle's mean position is determined using the mean Lagrangian velocity. Particle spreading (the displacement covariances) is expressed in terms of the state transition matrix that satisfies a time-varying dynamic equation whose coefficient matrix is the derivative of the mean Lagrangian velocity field. In the special cases of stationary velocity fields the transition matrix becomes the identical matrix, and our model reduces to the well-known model of Dagan <1984>. For nonstationary but unidirectional flow fields our model reduces to that of Butera and Tanda <1999> and Sun and Zhang <2000>. The validity of the transport model is examined by comparisons with Monte Carlo simulations for the following three cases: transport in a mean gravity-dominated flow, in an unsaturated flow with a water table boundary, and in a saturated-unsaturated flow. An excellent agreement is found between our model results and those from Monte Carlo simulations.