In this study, on the basis of the van Genuchten–Mualem constitutive relationship, we develop a general nonstationary stochastic model for transient, variably saturated flow in randomly heterogeneous media with the method of moment equations. We first derive partial differential equations governing the statistical moments of the flow quantities by perturbation expansions and then implement these equations under general conditions with the method of finite differences. The nonstationary stochastic model developed is applicable to the entire domain of bounded, multidimensional vadose zones or integrated unsaturated–saturated systems in the presence of random or deterministic recharge and sink–source and in the presence of multiscale, nonstationary medium features. We demonstrate the model with some two-dimensional examples of unsaturated and integrated unsaturated–saturated flows. The validity of the developed stochastic model is confirmed through high-resolution Monte Carlo simulations. We also investigate the relative contributions of the soil variabilities (KS, α, and n) as well as the variability in recharge Q to the pressure head variance. It is found that the pressure head variance is sensitive to these variabilities, in the order of n, α, KS, and Q. Though the variability of α and n is usually smaller than that of KS and Q, their effect on the pressure head variance should not be ignored.