Existing fractal studies dealing with subsurface heterogeneity treat the logarithm of the permeability K as the variable of concern. We treat K as a multifractal and investigate its scaling and fractality using measured horizontal K data from two locations in the United States. The first data set was from a shoreline sandstone near Coalinga, California, and the second was from an eolian sandstone <Goggin, 1988>. By applying spectral analyses and computing the scaling of moments of various orders (using the double trace moment method <Lavallee, 1991; Lavallee et al., 1992>), we found that K is multiscaling (i.e., scaling and multifractal). We also found that the so-called universal multifractal (UM) <Schertzer and Lovejoy, 1987> model (essentially a log-Levy multifractal), was able to reproduce the multiscaling behavior reasonably well. The UM model has three parameters: &agr;, &sgr;, and H, representing the multifractality index, the codimension of the mean field, and the distance to stationary multifractal, respectively. We found (&agr;=1.7, &sgr;=0.23, H=0.22) and (&agr;=1.6, &sgr;=0.11, H=0.075) for the shoreline and eolian data sets, respectively. The fact that &agr; values were less than 2 indicates that the underlying statistics are non-Gaussian. We generated stationary and nonstationary multifractals and illustrated the role of the UM parameters on simulated fields. Studies that treated Log K as the variable of concern have pointed out the necessity for large data records, especially when the underlying distribution is Levy-stable. Our investigation revealed that even larger data records are required when treating K as a multifractal, because Log K is less intermittent (or irregular) than K. ¿ 2000 American Geophysical Union