The global semivariogram and covariance of ln(K) can be modeled with a hierarchy of correlation structures corresponding to the organization of bedding within a sedimentary sequence. Such a model accounts for the spatial correlation of ln(K) within and across bedding units defined at one level. This is related to correlation of ln(K) at higher levels (larger scales) through the spatial correlation of indicator variables representing the proportions, geometry and juxtaposition patterns of the units at each 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 illustrated with real data from a site with glaciofluvial sand and gravel deposits.