Geostatistical theory has shown promise in dealing with issues of stability, uniqueness, and identity of estimates inherent in inverse problems of subsurface flow. Here the geostatistical method is extended to three-dimensional, unsteady flow in variably saturated porous geological media (the vadose zone) that are modeled using the Richards equation and the van Genuchten-Mualem constitutive equations. The saturated hydraulic conductivity, &agr;, and n parameters of this relationship are treated as spatially correlated, statistically independent, stochastic processes for representing heterogeneity of porous media. For given covariance functions of the parameters the adjoint-state sensitivity method is used to calculate first-order approximations for covariances of capillary pressure and moisture content and cross covariances between capillary pressure, moisture content, and the hydraulic properties. These covariances and cross covariances are then used in a successive linear estimator (SLE) to estimate the conditional means of the heterogeneous hydraulic property fields based on measurements of pressure and moisture content data. A sequential conditioning approach for our SLE was also applied to data sets collected at different sampling times during a transient infiltration event. This approach has the benefit of reducing the size of the matrices and so helps avoid numerical stability problems. On the basis of our study, pressure and moisture content data sets collected at later times of an infiltration event or during steady state flow were found to provide better estimates (smaller mean-square error compared to the true field) of the hydrological parameters of the vadose zone than data from very early times. ¿ 2000 American Geophysical Union