FlowPaths Model and Simulations

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.本文提出一种名为FlowPaths的数值反演方法，用于基于观测得到的速度场求解多孔介质的水力传导率。该速度场可通过多种途径推导获得，例如含水层中的电阻率测量结果，或实验室实验中的粒子图像测速(Particle Image Velocimetry, PIV)技术。该反演方法假设离散化的块中心顶点网格内存在稳态二维流动。为了从速度场求解每个网格块的水力传导率K(x,y)，研究人员开发了一种图论方法，以找到贯穿多孔介质的线性独立流动路径集合。每条流动路径均可生成一条水头差方程，连接高水头边界上的源顶点与低水头边界上可抵达的汇顶点；随后借助达西定律(Darcy's Law)，将水头差方程中的各项转换为水力传导率的等效表达式。随后可直接求解该方程组以得到K(x,y)。将反演结果与用于生成速度场的合成非均质地力传导率场进行对比后，结果表明该反演方法在多种场景下均具备高精度，最大相对误差小于百万分之一；同时该方法还被证明具备递推稳定性。上述结果证实，FlowPaths可应用于现场与实验室场景，从已知速度场中反演得到水力传导率参数。 已知的合成非均质地力传导率场均通过MATLAB生成。准二维流动矩阵的网格尺寸范围为4×4至16×16，每个单元格均为边长1cm的正方形。每个流动矩阵的非均质性由基于储层指数Vdp的高斯随机数生成器控制，Vdp以1.0×10^(-8) cm²的绝对渗透率为中心，取值范围为0.01至0.98，据此可生成的K(x,y)非均质性范围约为1.6×10^(-7) cm/s至最大流动矩阵对应的200 cm/s以上。比流量通过MATLAB中的交替方向隐式(Alternating Direction Implicit, ADI)格式计算得到，流动矩阵的对边设置为ΔH=1cm的恒定水头边界，另外两侧设置为无流动边界。每次运算的收敛容差为1×10^(-9) cm水头，通过对比水平计算与垂直计算的结果实现收敛判断。 单元间流动方向(q)被用于构建有向图，边的权重等于比流量q。利用最小权重有向生成树算法，可识别出各对向边界单元格（顶点）之间的有向边路径。随后将这些边路径转换为顶点路径，并定义新变量q*为相邻边权重的平均值，以撰写各路径上顶点间的水头差方程dH=Σdh。将各独立水头差方程组装为稀疏矩阵A，总水头差dH记为向量b，随后求解矩阵方程Ax=b以得到x，其中x为传导率K的倒数。将求解结果与初始数据进行对比以验证精度，此外还可通过一项不依赖初始数据的递推过程，额外验证算法的稳定性。 本文中使用了一组输入数据作为示例：一个4×4的网格，其K的粗略量级(Rough Order of Magnitude, ROM)范围为0.82 cm/s，以10.2 cm/day为中心。反演求解结果与递推稳定性检查结果均已公开。 另一组涵盖所有网格尺寸（4×4至16×16）的大规模输入数据集已单独构建完成。所有反演结果与初始误差检查数据均已包含在本数据集中。