In this study we develop a stochastic model for transient unsaturated-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 stochastic model developed is applicable to the entire domain of a bounded, multidimensional unsaturated-saturated system in the presence of random or deterministic recharge and sink/source and in the presence of multiscale, nonstationary medium features. The unsaturated and saturated zones are coupled through the water table, whose position is random in randomly heterogeneous porous media. The presence of the water table renders the flow moments strongly nonstationary even in the absence of medium nonstationary features. This finding is confirmed with Monte Carlo simulations. The resulting first two moments of flow may be used to approximate confidence intervals for the flow quantities. In addition, the integrated stochastic flow model provides the prerequisite quantities for realistically studying contaminant transport in unsaturated-saturated systems with a stochastic approach.