In this paper we present a complex network model based on a heterogeneous preferential attachment scheme to quantify the structure of porous soils. Under this perspective pores are represented by nodes and the space for the flow of fluids between them is represented by links. Pore properties such as position and size are described by fixed states in a metric space, while an affinity function is introduced to bias the attachment probabilities of links according to these properties. We perform an analytical study of the degree distributions in the soil model and show that under reasonable conditions all the model variants yield a multiscaling behavior in the connectivity degrees, leaving a empirically testable signature of heterogeneity in the topology of pore networks. We also show that the power-law scaling in the degree distribution is a robust trait of the soil model and analyze the influence of the parameters on the scaling exponents. We perform a numerical analysis of the soil model for a combination of parameters corresponding to empirical samples with different properties, and show that the simulation results exhibit a good agreement with the analytical predictions