An iterative co-conditional Monte Carlo simulation (IMCS) approach is developed. This approach derives co-conditional means and variances of transmissivity (T), head (φ), and Darcy's velocity (q), based on sparse measurements of T and φ in heterogeneous, confined aquifers under steady-state conditions. It employs the classical co-conditional Monte Carlo simulation technique (MCS) and a successive linear estimator that takes advantage of our prior knowledge of the covariances of T and φ and their cross-covariance. In each co-conditional simulation, a linear estimate of T is improved by solving the governing steady-state flow equation, and by updating residual covariance functions iteratively. These residual covariance functions consist of the covariance of T and φ and the cross-covariance function between T and φ. As a result, the non-linear relationship between T and φ is incorporated in the co-conditional realizations of T and φ. Once the T and φ fields are generated, a corresponding velocity field is also calculated. The average of the co-conditioned realizations of T, φ, and q yields the co-conditional mean fields. In turn, the co-conditional variances of T, φ, and q, which measure the reduction in uncertainty due to measurements of T and φ, are derived. Results of our numerical experiments show that the co-conditional means from IMCS for T and φ fields have smaller mean square errors (MSE) than those from a non-iterative Monte Carlo simulation (NIMCS). Finally, the co-conditional mean fields from IMCS are compared with the co-conditional effective fields from a direct approach developed by Yeh et al. [Water Resources Research, 32(1), 85-92, 1996].