A numerical approach for modeling unsaturated flow is developed for heterogeneous simulations of fractured tuff generated using a geostatistical method. Cross correlations of hydrologic properties and upscaling of moisture retention curves is discussed. The approach is demonstrated for a study of infiltration at Yucca Mountain