Most theoretical work on analyzing plume spreading at the field scale in a partially saturated heterogeneous formation assumes weak stationarity of velocity field. While this assumption is not applicable to the case of a bounded flow domain, the nonstationarity in fluctuations of unsaturated velocity fields is induced through the presence of the boundary conditions in the flow domain. In this work, we attempt to quantify the large-time macrodispersion in nonstationary unsaturated velocity fields caused by the presence of a fixed head boundary condition. Application of the perturbation-based nonstationary spectral approach leads to an analytical expression for the macrodispersion describing the field-scale dispersive solute flux in terms of the statistical moments of two formation parameters: the Gardener's parameter (α) and the saturated hydraulic conductivity (Ks). The results predicted from the expression indicate that the enhanced unsaturated plume spreading can arise from a larger correlation scale of ln Ks or ln α fields. In addition, the α-parameter takes the role in reducing the field-scale plume spreading.