A numerical inverse method, called FlowPaths, is presented to solve for the hydraulic conductivity of a porous medium using an observed velocity field that can be inferred, for example, from electrical resistivity measurements in aquifers or from particle image velocimetry (PIV) in laboratory experiments. The inverse method assumes steady, two-dimensional flow in a discretized grid of block-centered vertices. To determine the hydraulic conductivity, K(x,y), of each grid block from the velocity field, a graph-theoretical approach was developed to find a linearly independent set of flow paths through the porous medium. Each path generates a head-drop equation connecting a source vertex on the high head boundary to a reachable sink vertex on the low head boundary, and then each of the terms in the head-drop equation is converted to an equivalent expression for hydraulic conductivity using Darcy’s law. The system of equations is then solved directly for K(x,y). When results are compared to the synthetic heterogeneous hydraulic conductivity field used to generate the velocity field, the inverse method is demonstrated to be accurate, with a maximum relative error of less than one part per million, in a wide range of scenarios. The inverse method is also shown to be recursively stable. These results show that FlowPaths can be used in field and laboratory applications to find the hydraulic conductivity parameter from a known velocity field.
The known synthetic heterogeneous hydraulic conductivity fields are created using MATLAB. Grid sizes of the quasi-2D flow matrices range from 4x4 to 16x16, with each cell having a square dimension of 1cm. The heterogeneity of each flow matrix is controlled by a Gaussian random number generator based on a reservoir index number (Vdp, as a function of absolute permeability centered on 1.0x10^(-8) cm^2 ) that ranges from 0.01 to 0.98, producing ranges of K(x,y) heterogeneity from approximately 1.6x10^(-07) to over 200 cm/sec for the largest flow matrices. The specific discharges were computed using an ADI scheme in MATLAB, with constant head boundaries of delta H=1cm on opposing sides and no-flow boundaries on the other two sides of the flow matrices. The tolerance for completion of each run was 1x10^(-09)cm of head, comparing the horizontal computation to the vertical computation.
The directions of intercellular flow (q) are used to create directed graphs, with edges having weights equal to q. Directed edge paths between each of the opposing boundary cells (vertices) are identified using minimum weight directed spanning trees. These edge paths are then used to formulate vertex paths where a new variable, q*, is identified as the average of adjacent edge weights to write head drop equations, dH=sum(dh) between each vertex, over each path. The individual head-drop equations are assembled as a sparse matrix A, the total head drop dH=b. The matrix equation Ax=b is solved for x, where x is the reciprocal of the conductivity K. The results are checked against the initial data for accuracy, and an additional check is available for concluding whether the algorithm is stable by a recursive process that does not depend on the original data.
One set of input data was used as an example in the paper: a 4x4 grid with a range of rough order of magnitude (ROM) of K=0.82 cm/sec, centered at 10.2 cm/day. The results of the inverse solution and the recursive stability check are also available.
A large set of input data was separately created using all sizes of grids (4x4 to 16x16). All of the inverse results and the initial error checks are available in the data set.