We introduce a stochastic model of flow through highly heterogeneous, composite porous media that greatly improves estimates of pressure head statistics. Composite porous media consist of disjoint blocks of permeable materials, each block comprising a single material type. Within a composite medium, hydraulic conductivity can be represented through a pair of random processes: (1) a boundary process that determines block arrangement and extent and (2) a stationary process that defines conductivity within a given block. We obtain second-order statistics for hydraulic conductivity in the composite model and then contrast them with statistics obtained from a standard univariate model that ignores the boundary process and treats a composite medium as if it were statistically homogeneous. Next, we develop perturbation expansions for the first two moments of head and contrast them with expansions based on the homogeneous approximation. In most cases the bivariate model leads to much sharper perturbation approximations than does the usual model of flow through an undifferentiated material when both are applied to highly heterogeneous media. We make this statement precise. We illustrate the composite model with examples of one-dimensional flows which are interesting in their own right and which allow us to compare the accuracy of perturbation approximations of head statistics to exact analytical solutions. We also show the boundary process of our bivariate model is equivalent to the indicator functions often used to represent composite media in Monte Carlo simulations.