In this study we develop a first-order, nonstationary stochastic model for steady state, unsaturated flow in randomly heterogeneous media. The model is applicable to the entire domain of a bounded vadose zone, unlike most of the existing stochastic models. Because of its nonstationarity, we solve it by the numerical technique of finite differences, which renders the flexibility in handling different boundary conditions, input covariance structures, and soil constitutive relationships. We illustrate the model results in one and two dimensions for soils described by the Brooks-Corey constitutive model. It is found that the flow quantities such as suction head, effective water content, unsaturated hydraulic conductivity, and velocity are nonstationary near the water table and approach stationarity as the vertical distance from the water table increases. The stationary limits and the critical vertical distance at which stationarity is attained depend on soil types and recharge rates. The smaller the recharge rate is, the larger the critical distance; and the coarser the soil texture is, the smaller the distance. One important implication of this is that the existing simpler, gravity-dominated flow models may provide good approximations for flow in vadose zones of large thickness and/or coarse-textured soils although they may not be valid for vadose zones of fine-textured soils with a shallow water table. It is also found that the vertical extent of a domain where nonstationarity is important may be estimated by solving the one-dimensional Richards equation for mean head with average soil properties and appropriate boundary conditions. On the basis of the mean head, one may then determine whether the full, nonstationary model must be solved or whether a simpler, gavity-dominated model will suffice. The flow quantities are also nonstationary in the horizontal direction near the lateral boundaries, as found for flow in saturated zones.