Stochastic imaging provides geologically realistic solutions to aquifer uncertainty quantification problems. Higher-order statistics for representing complex priors in stochastic-imaging problems are typically borrowed from training images, which may bias outcomes if the training images are unrepresentative of the desired structures. A Markov random field (MRF)-based stochastic-imaging algorithm, the spatial statistics of which are adaptively inferred directly from hydrogeophysical measurements without training images, is presented. For convenience of implementation, MRF models are commonly specified with pairwise interactions, limiting their ability to model complex architectures. MRF model specifications involving higher-order interactions are considered here. In the proposed algorithm, the lithologic structure of the aquifer and the hydraulic conductivities within the identified lithologies are simultaneously estimated in conformance to data-unconditioned and data-conditioned higher-order statistics while honoring the hydrogeophysical data sets. High reconstruction accuracy rates with nonsmooth, geologically realistic aquifer heterogeneities are reported in a hypothetical trinary hydrofacies aquifer characterization problem with various combinations of concentration and electrical-resistivity conditionings, illustrating the value of the different data types.