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