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: alpha, sigma, and H, representing the multifractality index, the codimension of the mean field, and the "distance" to stationary multifractal, respectively. We found (alpha = 1.7, sigma = 0.23, H = 0.22) and (alpha = 1.6, sigma = 0.11, H = 0.075) for the shoreline and eolian data sets, respectively. The fact that alpha 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.