SimPEG+Python：地下水流有限体积建模与稀疏求解实践

简介：这是一套基于SimPEG框架编写的地下水流方程模拟源码包，面向从事地球物理、水文地质以及环境工程研究的Python开发者，可用于学习地下水流正演与参数反演的完整建模流程。压缩包共包含20个文件，其中9个Python脚本、2个RST文档、Makefile及若干配置和说明文件，整体仅45KB，结构非常紧凑；同时附有独立文档与多个示例目录，便于对照学习与二次开发。已有492人浏览学习，这一项目特别适合需要快速掌握SimPEG.FLOW的开发者：从中可了解达西定律驱动的饱和/非饱和水流模拟、多种边界条件设置，以及基于观测数据的参数反演思路。借助源码中的网格构建、优化算法和可视化代码，读者能直接复用这些模块来处理污染物扩散、水源保护等实际问题，提升地下水建模效率。

1. 用 SimPEG 框架编写地下水流方程：先别指望“开箱即用”

当一份资源标题写着“用 SimPEG 框架编写地下水流方程”时，直觉上会觉得它又是一套开箱即用的地下水模拟器。实际不是。SimPEG 只给你网格、梯度算子、散度算子和稀疏求解所需的基础工具，达西定律、水头方程、边界条件全部要用 Python 和 NumPy 逐行组装。这种“半框架”的好处是：你能看见方程每个系数落在矩阵哪一区，参数怎么改心里有数。simpegflow 这份资源把一维渗流、二维抽水井流场和参数化封装成可直接跑的脚本，适合想跨进地下水流模拟、又不想被重型软件黑匣子劝退的 Python 从业者。新手能照步骤敲出结果，熟手能在源码基础上改边界和分层。

2. 把地下水流方程“摆”上网格：达西定律、有限体积与边界条件的对应关系

2.1 控制方程与符号约定：水头才是主角

承压含水层中的稳态地下水运动，由达西定律描述：单位面积通量 q = -K ∇h。这里的 h 是水头，单位 m；K 是水力传导系数，单位 m/d；q 是达西通量，单位 m/d。把 q 代入质量守恒方程，就得到稳态地下水流方程：

∇·(K∇h) = 0

如果有抽水井或入渗补给，右端加一个源汇项 Q。写代码之前，我习惯先确定一个符号约定：K定义在哪里，h定义在哪里，Q的正负号怎么算。simpegflow 里默认 h 放在单元中心，K 也放在单元中心，面通量 q 在面上参与计算。这套约定和 SimPEG 的离散算子天然匹配，下面会讲到为什么不能随手换。

注意：simpegflow 默认求解稳态方程。瞬态问题要把储水系数 S 和单位体积质量矩阵组装进去，这部分在资源里只给了注释，没展开主流程，我建议新手先把稳态吃透再碰瞬态。

2.2 有限体积与算子的对应：梯度在面上，散度在单元内

SimPEG 框架的底层是有限体积，核心对象是张成网格（TensorMesh）。它与常规有限差分的关键区别在于：场变量被离散到“单元中心”，而导数算子和材料属性被离散到“面中心”。因此在组装过程中，你要区分两个完全不同的几何对象：

• 单元中心（cell center）：放水头 h、水力传导系数 K 这类“块状”物理量；
• 面中心（face center）：放通量、梯度、法向流速这类“通过面交换”的量。

这张对应关系决定了 SimPEG 中两个算子的形状。cell_gradient_cc接收单元中心的值，输出面中心的梯度；face_divergence接收面中心的通量，输出单元中心的散度。二者组合起来就得到水头到“单元内水量净收支”的映射：

A = Div * diag(K_face) * Grad

其中K_face是把单元中心 K 插值到面中心之后的对角阵。这个 A 矩阵的每一行，对应一个单元的水量平衡方程。把源汇项放到右侧，就得到线性系统：

A h = b

这里有个习惯问题：SimPEG 的张成网格默认边界是零通量（封闭边界），所以 Neumann 边界你不用额外处理，不做任何修改就是绝源边界。这一点第一次上手容易漏，我见过不少人在边界上手动清零，反而把默认的封闭性破坏了。

2.3 边界条件的三种矩阵化写法

边界条件在有限体积里不是“设置参数”，而是“改写矩阵条目”。simpegflow 用到的三类边界：

写代码时，最容易出错的是第二类和第三类的叠加。边界单元同时被定通量和定水头约束时，矩阵行只能保留一种主约束；通常的做法是保留 Cauchy 或 Dirichlet，而把通量约束写进残差项。这个取舍在二维模型里会直接影响井点附近的收敛，后面避坑章再展开。

3. 一维稳态渗流实战：30 行代码把分层含水层水头剖面算出来

3.1 两层介质、左右定水头的模型设计

先跑通一维模型，是验证框架和算子理解最便宜的办法。模型设定如下：一维含水层总长 1000 m，左段 0 到 400 m 是砂质层，水力传导系数 K1 = 10 m/d；右段 400 到 1000 m 是粘土质层，K2 = 0.5 m/d。左右边界都是定水头：左侧 25 m，右侧 10 m。水从左向右流动。

我不想把这个案例做成“背代码”，而想把建模过程拆成三步：建网格、赋属性、组装求解。网格划分可以选择均匀分 100 个单元，每段 10 m。尽管两层介质交界面在 400 m 处正好落在单元边界上，离散化时不会出现跨层单元，这是刻意为之，方便检查结果。

3.2 从网格到矩阵：完整的 Python 实现

simpegflow 中的一维脚本核心代码如下。注意K_cc是单元中心的属性向量，要用aveCC2F插值到面中心，再参与组装。

import numpy as np import scipy.sparse as sp from scipy.sparse.linalg import spsolve from discretize import TensorMesh # 1. 建网格：100 个均匀单元，每个单元 10 m mesh = TensorMesh([np.ones(100) * 10.0]) n_cells = mesh.n_cells # 2. 定义单元中心水力传导系数：前 40 个单元是砂层，后面是粘性层 K_cc = np.ones(n_cells) * 0.5 # 先赋粘性层 K_cc[:40] = 10.0 # 0~400 m 砂层 # 3. 核心操作：单元中心属性插值到面中心 K_face = mesh.aveCC2F * K_cc # 4. 组装散度-通量矩阵 Grad = mesh.cell_gradient_cc # (n_faces, n_cells) Div = mesh.face_divergence # (n_cells, n_faces) A = Div @ sp.diags(K_face) @ Grad # 5. Dirichlet 边界处理：替换首末两行 A = A.tolil() b = np.zeros(n_cells) A[0, :] = 0 A[0, 0] = 1 b[0] = 25.0 A[-1, :] = 0 A[-1, -1] = 1 b[-1] = 10.0 # 6. 稀疏求解 h = spsolve(A.tocsr(), b)

这段代码里，mesh.aveCC2F是单元中心到场面的平均算子，它处理的是“从块中心到面中心”的插值。cell_gradient_cc返回面中心处的梯度值，维度和面数一致。face_divergence再把面通量聚合到单元中心。这三个算子本身是 SimPEG 框架给好的，我们要做的只是用矩阵乘法把它们组合成 A。组装顺序是：先插值 K，再乘梯度，最后散度；不能换乘顺序，否则稀疏矩阵作用对象就反了。

行替换是处理 Dirichlet 边界最省事的方式，但会在边界单元上破坏原来的通量守恒方程。严格做法是把边界通量对应的项移到右端，等避坑章会讲这事的定量影响。

3.3 结果校验：通量连续性和残差检查

跑完spsolve后不能立刻信结果，至少要做一个“通量连续性检查”。在地下水流里，稳态条件下流量在任意断面上都应相等。我通常用这四行代码验证：

# 由水头结果反算面通量 q_face = -K_face * (Grad @ h) # 取出左右边界面上的通量，再取交界面两侧单元的通量 q_left = -K_face[0] * (h[0] - 25.0) / 10.0 q_right = -K_face[-1] * (10.0 - h[-1]) / 10.0 print(q_left, q_right) # 计算内部单元残差：内点行应接近 0 res = (Div @ (K_face * (Grad @ h)))[1:-2] print(np.abs(res).max())

这个检查能暴露两类问题：如果q_left和q_right明显不相等，说明边界替换影响了内部解；如果内部残差不是极小值，说明算子组合的方向或符号写反了。常见做法里，残差考察的是“面通量散度回到单元里”是否为零，而不是看水头本身。simpegflow 把这段检查直接写在脚本末尾，我建议保留，后续所有二维模型都值得复用。

4. 二维抽水井流场：网格尺寸、导水系数分层与后处理

4.1 模型设定与网格离散：区域、井位、边界

二维模型才能体现地下水流模拟的实用性。设计一个 2000 m × 1000 m 的矩形承压含水层，x 方向 0 到 2000 m，y 方向 0 到 1000 m。上层（y > 400 m）水力传导系数 K1 = 20 m/d，下层 K2 = 5 m/d，模拟砂砾层覆盖在粉砂层上的情况。在区域中心 (1000, 500) 处设一口抽水井，流量固定 1200 m³/d，即源项为负。矩形四周边界都为定水头 30 m，相当于一个被湖泊包围的含水层。

网格划分我一般不在全区域用均匀边长，而是让中心附近加密、外边界方向逐渐放大。二维TensorMesh需要分别给定 x 和 y 方向的距离数组，例如：

hx = [(50.0, 20, -1.3), (10.0, 20), (50.0, 20, 1.3)] hy = [(50.0, 10, -1.3), (10.0, 10), (50.0, 10, 1.3)] mesh = TensorMesh([hx, hy])

三元素的元组写法是“从边界向内部加密”的常用做法：第一项是起始单元尺寸，第二项是单元个数，第三项可选，正数为向外边长放大因子，负数为向内。这样在井点附近分辨率高，外围单元数量少，节省稀疏矩阵内存。

4.2 二维组装与求解：边界行替换的高效写法

二维情况下边界单元数量变多了，逐行tolil逐个替换不现实。更高效的做法是先用布尔索引找出所有贴边单元，一次批量替换矩阵行。核心代码如下：

import numpy as np import scipy.sparse as sp from scipy.sparse.linalg import spsolve from discretize import TensorMesh # 网格生成：中心加密，边界放大 hx = [(100.0, 15, -1.3), (20.0, 20), (100.0, 15, 1.3)] hy = [(100.0, 10, -1.3), (20.0, 10), (100.0, 10, 1.3)] mesh = TensorMesh([hx, hy]) n_cells = mesh.n_cells cc = mesh.cell_centers # 分层导水系数：y 大于 400 m 为高渗层 K_cc = np.ones(n_cells) * 5.0 K_cc[cc[:, 1] > 400.0] = 20.0 # 插值到面中心并组装 K_face = mesh.aveCC2F @ K_cc Grad = mesh.cell_gradient_cc Div = mesh.face_divergence A = Div @ sp.diags(K_face) @ Grad # 源汇项：抽水井落在最接近中心的单元 well_index = np.argmin((cc[:, 0] - 1000.0) ** 2 + (cc[:, 1] - 500.0) ** 2) cell_vol = mesh.cell_volumes[well_index] b = np.zeros(n_cells) b[well_index] = -1200.0 / cell_vol # 四周边界水头 30m bc_cells = ( np.isclose(cc[:, 0], 0.0, atol=50.0) | np.isclose(cc[:, 0], 2000.0, atol=50.0) | np.isclose(cc[:, 1], 0.0, atol=50.0) | np.isclose(cc[:, 1], 1000.0, atol=50.0) ) A = A.tolil() for i in np.where(bc_cells)[0]: A[i, :] = 0.0 A[i, i] = 1.0 A = A.tocsr() b[bc_cells] = 30.0 h = spsolve(A, b)

这里的.atol=50.0是根据边界附近最大单元尺寸选的容差。因为网格在边界区域用了放大系数，最外层单元中心到边界的距离可能接近几十米，np.isclose容差不足会漏掉边界单元。参数-1200.0 / cell_vol是源汇项的单位换算：井流量 1200 m³/d，除以单元体积（单位 m³，二维视为 1m 厚度）之后，才能放进稳态方程的右端。

组装完成后，后处理建议用pcolor画水头等值面，再把网格边界叠加上去。重点观察井点周围的水头降落漏斗，以及分层界面处等值线间距的变化。高渗层水头变化平缓，低渗层梯度增大，这是判断分层赋值有没有生效的直观特征。

4.3 参数敏感性：K 倍率与水头降深的量化关系

只画出一张图，还不能证明模型可靠。我会固定井流量和边界条件，把上下层 K 同时放大 2 倍、5 倍，观察井点水头的响应。下图分三组运行后记录井点水头下降值：

井点水头随 K 放大而上升，符合“导水能力增强、水头损失减小”的定性预期。但要注意，水头变化幅度不是线性的，说明抽水井附近的局部几何效应参与了调节。这种量化表比单张云图有说服力，也便于你判断网格离散化是否已经收敛。

5. 避坑：SimPEG 地下水流建模最容易翻车的 5 处细节

5.1 现象一：把单元中心的 K 直接当作面通量，解出来水头剖面呈锯齿状

• 现象：一维模型求解后，水头曲线不是光滑斜线，而在每个单元之间上下波动，看起来像“锯齿”。
• 原因：直接把K_cc放进sp.diags，没有先用aveCC2F插值。有限体积的通量计算发生在面中心，而单元中心的 K 不能直接用于面梯度；未经插值会引入半个单元偏移，导致一阶误差变成符号交替的振荡。
• 解决：强制使用mesh.aveCC2F @ K_cc得到K_face。验证方式是打印K_face.shape，确认它的长度等于面数而不是单元数。任何形状对不上，都是属性位置放错了。

5.2 现象二：Dirichlet 只改矩阵行，右侧未写边界水头，结果整场偏低

• 现象：边界处水头设置成 30 m，求解完内部水位全部比预期低几米，边界单元正确但相邻单元错位。
• 原因：替换矩阵行只约束了未知量 h，但右侧向量 b 在对应行没有写入边界水头值。稀疏线性系统求解的是 A h = b，如果 b 是零，那么这一行变成 0 = 0 或 h = 0，取决于你的替换方式。
• 解决：批量替换时，先把b[bc_cells] = 30.0写进代码，再执行行替换；或者更安全的写法是先用全零向量初始化 b，再在替换循环里同步设置b[i] = 30.0。我在写二维代码时就吃过这个亏，现在每提交一次运行都检查np.unique(h[bc_cells])。

5.3 现象三：交界面的面插值把两层 K 平均出“平滑值”，分层边界失效

• 现象：两层导水系数差 40 倍，但模拟出的水头等值线在界面附近过渡得很平滑，看不到明显折射。
• 原因：aveCC2F在交界面两侧会做算术平均，把一个高渗单元和低渗单元的 K 平均成中间值，这不符合层状介质“边界两侧物理性质突变”的实际情况。
• 解决：要么在建模时把网格线刻意对准分层界面，让交界面不跨过单元；要么给aveCC2F换成调和平均方案（比如先算术平均再倒数平均的谐波平均），但 SimPEG 默认算子不带这个功能，需要自己写一个权重矩阵。simpegflow 的二维例子选的是前一种方案，靠网格线对齐。

5.4 现象四：命令能跑通、结果却整体偏大，问题往往出在单位制

• 现象：所有步骤都和示例一致，但水头数值偏大上百倍，井点降落漏斗不自然。
• 原因：K 用的是 m/d，但井流量用了 m³/s，或单元尺寸用了 km。有限体积的所有几何量都进入矩阵，单位不一致会让通量项和源汇项差好几个数量级。
• 解决：建模前固定一套单位制，建议全部用“米 + 天”组合：K 单位 m/d、井流量单位 m³/d、坐标单位 m。写代码时把单位写在参数名注释里，比如K_cc = 20.0 # m/d。这份资源的一维脚本开头有一行单位约定注释，别删。

5.5 现象五：几十万单元的二维模型直接 spsolve，内存当场爆炸

• 现象：加密网格到 20 万单元以上，spsolve长时间无响应，内存占用超过 16 GB。
• 原因：A.tolil()后的行替换会让稀疏结构变复杂，叠加spsolve的稠密分解中间量，内存峰值远比 A 本身大。
• 解决：边界替换改成csr行索引切片赋值，不要tolil后逐行改；对大模型先尝试scipy.sparse.linalg.cg，并给rtol一个宽松的初始值，比如 1e-6。如果网格可以粗化，先用 2 万单元跑通，再加密到 20 万验证结果差异小于 0.5%，这是最简单也最有效的后悔药。

6. 验算技巧：把解析解固化成回归测试，网格怎么加密都不心虚

二维抽水模型没有现成的解析解可套，但一维均匀介质有。最经典的验证案例是：单层均匀含水层，长 L，左右定水头 h0 和 h1，稳态水头分布理论解是直线 h(x) = h0 + (h1 - h0) * x / L。这个解简单到几乎“不像验证”，但却是检验整个算子链是否纯净最灵敏的探针。

我习惯在 simpegflow 的一维脚本里加一段自动回归测试：

L = 1000.0 x_cc = mesh.cell_centers h_analytic = 25.0 + (10.0 - 25.0) * x_cc / L h_numeric = solve_1d(mesh, K_cc, 25.0, 10.0) max_err = np.abs(h_numeric - h_analytic).max() assert max_err < 1e-6, f"数值解偏离解析解 {max_err:.2e}"

这段时间发现的教训是：矩阵组装里的任何一个小偏差，都会让数值解偏离直线，尤其是插入aveCC2F和 Dirichlet 行替换这两个环节。我第一次跑通二维模型之前，先拿这个解析解反复校验网格、算子、边界三类代码，确认max_err在 1e-6 以下，才敢去调抽水井参数。从那以后，我每次拿到新的网格设置或新的源汇项格式，都强制走一遍这个解析解回归，再转二维；这个习惯帮我挡掉了至少三次边界替换的隐性 bug。二维模型还可以另加一个对称校验：井位在中心时，水头场关于过井点的 x、y 对称线应是对称的，用np.allclose(h_reshaped, h_reshaped.T)就能快速检查。希望这些验算步骤帮你在 SimPEG 框架下少走一点弯路。

本文还有配套的精品资源，点击获取