The spatial covariance of ln(k) can be modeled with a hierarchy of covariance structures corresponding to the organization of bedding within and among the lithofacies of a sedimentary sequence. Such a model accounts for the spatial correlation of ln(k) within and across bedding units defined at any one level (Ritzi et al., 2004). This is related to correlation of ln(k) at a higher level (larger scale) through the spatial correlation of indicator variables representing the proportions, geometry, and juxtaposition patterns of the units at the lower level. In this paper the fitting of the components of the hierarchical model, written as nested functions, is considered in developing a hierarchical covariance model for use in estimation, simulation, or analytical derivation of macrodispersivity models. The least-squares criterion, along with parameter prior information and other weighted constraints, is used as the objective function of the inverse problem, which is solved by the Gauss–Newton–Levenberg–Marquardt method. The method is tested on synthetic data and illustrated with real data from a site with glaciofluvial sand and gravel deposits.