Most existing stochastic models are developed for unimodal porous media that may be well characterized with only the first two statistical moments. However, the distribution of hydraulic properties, such as hydraulic conductivity, may possess a multiplicity of modes; thus the first two moments may not be adequate to characterize properties of such porous media. In turn, the stochastic models developed for unimodal porous media may not be applicable to flow and transport in a multimodal heterogeneous porous medium. This study investigates under what circumstances the second-order moment-based stochastic models are applicable to multimodal heterogeneous porous media. We assume that a porous medium is composed of a number of materials (categories), each of which may have a different mean, variance, and correlation scale. The distribution of materials in the domain is characterized by indicator random variables. We first derive analytical expressions for the mean and covariance of the log saturated hydraulic conductivity (ln Ks) of the multimodal porous medium in terms of categorical proportions, transition probability among categories, and covariances of indicator random variables. We express the covariance in terms of the statistics of materials in the porous medium, which allows us to accurately evaluate the variance and the correlation length of the composite ln Ks field. We then solve the second-order moment equations for the equivalent unimodal field with an exponential covariance with a single correlation scale computed for the composite field. On the other hand, we conduct two sets of Monte Carlo simulations: one with multimodal random fields, and the other with equivalent unimodal random fields. Examples for porous media with two materials are given. Numerical experiments show that a bimodal ln Ks field may be well approximated by an equivalent unimodal field when one of the two modes is dominant, under which condition the applicability of the second-order moment-based model is subject to the same limitation of relatively small variance as that for unimodal fields. When the bimodal distribution has two more or less equally important modes, although it cannot be adequately represented by an equivalent unimodal distribution, the second-order moment-based stochastic model seems to be applicable to systems with larger composite variances than it does for an one-mode-dominant distribution.