地下水和地热系统数值建模简介
质量、能量和溶质的基本原理
多孔弹性岩石中的运输
J. Bundschuh 和 M.C. Suárez Arriaga
多物理场建模
第 2 卷
地下水和地热系统的数值建模简介
封面
主图代表一个传导-对流地热系统，流体、热量和溶质通过孔隙和裂缝输送。简化的地质模型对应于墨西哥 Los Azufres 的地热田。
方程（从上到下）是以下各项的数学模型：
1） 压力、密度和渗透率方面的地下水流方程（方程 4.16）。
2） Terzaghi 有效应力和孔隙流体压力（方程 2.69b）。3） 多孔弹性应变、应力和多孔弹性系数的逆矩阵（方程 2.77b）。
4） 根据机械应力和热应力表示的孔隙流体压力（方程 2.120）。
5） 具有色散和对流的一般溶质输运方程（方程 4.111）。
6） 傅里叶热传导定律（方程 4.79）。
7） 具有传导和对流的一般热流方程（方程 4.104a）。
多物理场建模
系列编辑器
约亨·邦德舒
应用科学大学， 应用研究所， 德国卡尔斯鲁厄皇家理工学院 （KTH）， 瑞典斯德哥尔摩
马里奥·塞萨尔·苏亚雷斯·阿里亚加
米却肯大学 UMSNH 物理与数学科学学院应用数学与地球科学系，墨西哥米却肯州莫雷利亚
国际标准刊号：1877-0274
第 2 卷
地下水和地热系统数值建模简介
多孔弹性岩中的质量、能量和溶质传输基础
约亨·邦德舒
应用科学大学， 应用研究所， 德国卡尔斯鲁厄皇家理工学院 （KTH）， 瑞典斯德哥尔摩
马里奥·塞萨尔·苏亚雷斯·阿里亚加
米却肯大学 UMSNH 物理与数学科学学院应用数学与地球科学系，墨西哥米却肯州莫雷利亚
114
数值建模简介
Los Humeros 地热储层的一般特征：
平均储层压力 = 12.5 MPa 平均温度 = 310C 初始自然状态 = 压缩液体 最高温度 = 400C 估计电容量 = 75 MW 最小容积 = 19 km
3.8.3
Los Azufres， 米却肯州
1988 年对 Los Azufres 储层岩心的岩石物理特性进行了检查。从 18 个孔中收集了 24 个核心的数据。在 Az-19 井进行了一次扩散率和比热测量：在 250C 和 8.0 MPa 下，热扩散率为 0.66 × 10m /s，比热 = 1165 J/kg/C。测量是在 Terra Tek 岩石力学实验室进行的（Contreras 等人，1988 年）。岩石物理数据如表 3.8 所示。
Los Azufres 地热储层的一般特征：
平均储层压力 = 9.2 MPa （N）， 10.0 MPa （S） 平均温度 = 301C （N）， 311C （S） 初始自然状态 = 液体 （N），两相 （S） 最高温度 = 360C （S） 估计电容量 = 300 MW （N）， 200 MW （S） 最小体积 = 49 km （N）， 23 km （S） （N 和 S 代表 Los Azufres 地热田的北部和南部， 分别）。
表 3.8.
来自 Los Azufres 地热田的一些岩石物理参数（数据改编自 Contreras 等人，1988 年）。
Well
深度 （m）
ϕ (%)
ρ （千克/米）
k （mD）
千（W/m/C）
c（J/kg/C）
AZ-01 系列
1825.0–1829.6
2.6
2720
0.0010
–
–
AZ-03 系列
600.0–605.5
14.8
2300
0.0035
1.68
–
AZ-03 系列
1874.0–1880.0
13.2
2560
0.1773
1.84
–
AZ-03 系列
2117.0–2119.7
2.1
2740
0.0013
1.99
–
AZ-04 系列
1000.0–1000.5
12.6
2430
0.0018
1.56
–
AZ-05 系列
600.0–600.5
23.2
2080
0.0017
1.17
–
AZ-05 系列
1160.0–1165.0
11.9
2380
0.1513
–
–
AZ-08 系列
800.0–804.0
7.8
2590
0.1235
2.34
–
AZ-09 系列
1705.0–1710.0
2.6
2660
2.224
–
–
Az-10
1004.0–1005.0
4.7
2660
0.0013
–
–
Az-19
1000.0–1005.0
15.5
2290
0.015
1.97
1165
Az-20
650.0–654.5
13.1
2260
0.0018
1.58
–
Az-20
1600.0–1603.0
4.7
2660
0.0015
1.71
–
Az-22
800.0–805.0
9.9
2450
0.0017
2.17
–
Az-25
671.0–675.0
14.5
2300
0.0018
1.75
–
Az-26
596.0–601.0
2.6
2610
0.0020
2.20
–
Az-26
1002.0–1007.0
10.4
2410
0.4010
1.55
–
Az-29
400.0–402.0
20.1
2070
0.0410
1.05
–
Az-29
2496.0–2496.2
0.7
2810
–
–
–
AZ-41 系列
600.0–605.0
16.3
2360
0.0013
–
–
阿兹-46
802.0–805.0
7.4
2510
–
–
–
阿兹-47
2962.0–2964.0
2.1
2760
0.0020
1.89
–
Az-48
2678.5–2684.5
1.0
2840
–
–
–
AZ-50 系列
1133.0–1136.0
8.9
2470
0.0100
1.52
–
平均
–
9.3
2500
–
1.75
1165
第 4 章
流体流动、热量和溶质输送
“我第一次尝试解决瞬态问题是采用 Thiem 的受限条件方程，将其应用于具有自由表面的地下水体，并想象从储存中抽出的水奇迹般地被输送到 Thiem 锥体的外缘，并从那里渗入井中。然后，我可以计算这个 Thiem's 锥的体积，乘以我当时所说的比产率，将其等于泵送速率乘以时间，当然，还可以得到一个以时间表示的 Thiem's 锥外半径方程。然后，这可以替换为 Thiem 方程中的外半径值。你得到的是一个瞬态方程，它与现在的非平衡方程相同，只是 u 的井函数只包含对数项。常数和长幂级数都缺失了。
亨利·达西 （1803–1856）
4.1
流体的质量守恒
多孔体积中的总流体质量可以使用其密度分布来计算，就像在 2.1.3 节中所做的那样;在任何差分多孔体积 dV 中：
以 dV 为单位的流体质量：dM= ρdV= φ ρdV
⇒
M=
∫
V
φ ρ（r， t）dV
r = （x， y， z） 是流体粒子的位置矢量，单位为 V= φV;t ≥ 0.
(4.1)
该方程假设孔隙流体密度取决于空间和时间。如果流体的质量是守恒的，则它在运动中必须保持恒定。因此，前一个积分 （4.1） 的平流或总导数（方程 2.4）为零（方程 2.6a）：
D M
Dt
=
D
Dt
∫
V
φ ρ（r， t）dV =
∫
V
(
∂(ϕ ρ)
∂t
+ ∇ ·（φ ρv ）
)
dV = 0
(4.2)
该方程是多孔岩石中流体质量守恒原理的积分形式，适用于岩石连续体的任何部分。矢量 v = （v， v ， v ） 是流体粒子的场速。图 4.1 显示了具有体积模量 Vand 孔隙体积 V 的差异多孔岩石体积的示意图。体积 V 可以是多孔岩石的任意部分。例如，设 Vbe 一个微分体积 （REV） V ≈ dV，因此方程 （4.2） 的连续被积函数必须等于零：
∂(ϕ ρ)
∂t
+ ∇ ·（φ ρv ） =
D（φ ρ）
Dt
+ (ϕ ρ) ∇ ·v = 0
(4.3)
∗ 资料来源：Robert R. White 和 Alfred Clebsch：CV Theis：人及其对水文地质学的贡献。收录于：C.V. Theis 对地下水水文学的精选贡献，以及他的生活和工作回顾，美国地质调查局，供水论文 2415,1994 年。
115
116
数值建模简介
dV
v
谷物
V
V
dS
dV
v
固体
毛孔
dV
图 4.1.
简化多孔岩石的体积差异 dV 和表面 dS;v 是流体粒子的场速度。
公式 （4.3） 也称为流体的连续性方程，它以偏微分方程的形式表示相同的质量守恒原理。在适当的条件下，该定律可以与流体的体积变形直接相关（有关详细信息，请参见第 2.1.3 节）：
D（φ ρ）
Dt
= −(ϕ ρ) ∇ ·v
⇒
∇ ·v =
−1
ϕ ρ
D（φ ρ）
Dt
;
但 Vρ= M
⇒
−V
M
D
Dt
(
米
V
)
= −V
D
Dt
（V ）= V
(
V f
D V
Dt
)
=
1
V
(
D V
Dt
)
因此：
with： ε=
V
V
= ∇ ·u ,
v =
∂ u
∂t
⇒
∇ ·v =
1
V
(
D V
Dt
)
=
D ε
Dt
(4.4)
因此，流体速度的散度等于流体体积 V 的变化率。变量 ε 定义为液相的体积变形（膨胀或压实）。如果流体和孔隙率不随时间变化：
∂(ϕ ρ)
∂t
= 0
⇒
∇ ·（φ ρv ） = 0
(4.5)
如果流体是不可压缩的并且孔隙率是恒定的：
ϕ ρ∇ ·（v ） = 0
⇒
∇ ·v= 0
(4.6)
前面的所有方程都是同一流体质量守恒原理的不同形式。
最后一个公式 （4.6） 是连续性方程的最简单形式。如果在储液罐的一个或多个部分提取或注入 （q） 流体质量：
∫
V
(
∂(ϕ ρ)
∂t
+ ∇ ·（φ ρv ）
)
dV =
∫
V
q （r， t）dV
⇔
∂(ϕ ρ)
∂t
+ ∇ ·（φ ρv ） = q
[
kg
ms
]
(4.7)
流体流动、热量和溶质输送
117
其中 q表示从储层（汇或井：q < 0）中流体的体积抽取或在钻井的特定位置将流体注入储层（来源：q> 0）。
4.2
流体流动的一般模型：NAVIER-STOKES 方程
1822 年，伟大的法国数学家和工程师克劳德·路易斯·纳维耶 （Claude Louis Navier） 在流体流动方程中引入了粘度的影响。1845 年，英国数学家和物理学家乔治·加布里埃尔·斯托克斯 （George Gabriel Stokes） （图 4.2） 以类似于我们今天使用的形式推导出了相同的方程。从那时起，Navier-Stokes 方程就被用于解决科学和工程的许多分支中的流体流动问题。这些方程已被证明是“物理学中所有偏微分方程中最具挑战性的方程”之一（Brodkey 1967）。它们在流体内部包含摩擦效应，因此它们在物理上比 1822 年之前使用的欧拉方程更真实。
Navier-Stokes 方程描述了真实粘性流体的一般流动。他们起源
来自流体动力学基本定律和牛顿流体中的粘性力的结合。它们的构造方式如下。我们首先定义作用在流体上的主要函数和系数：σ 是应力张量，λ 和 μ 是流体粘度系数，p 是流体压力，Dis 是流体速度空间变化的对称张量。牛顿流体的行为由以下方程定义 （Germain 1973）：
张量形式：σ = −p I + λD I + 2 μD
组分形式：σ= （−p + λD ）δ+ 2 μD
D =
1
2
(
∂v
∂x
+
∂v
∂x
)
,
D=
∂v
∂x
+
∂v
∂y
+
∂v
∂z
= ∇ ·v
(4.8)
其中 λ 是膨胀（或压缩）的粘度，μ 是剪切流体动力学粘度。组成系数 K = λ+2μ/3 是流体的体积模量或体积动态粘度。这三个系数的单位是 [1 Pa · s = 1 kg/m/s = 1 g/cm/s = 1 泊]。流体动力学的基本定律是：
ρa = ∇ ·σ + F
(4.9)
其中向量 A 是流体粒子的加速度，F 是施加到流体上的外部矢量力（参见第 2.2.10 节）。计算等式 （4.8） 中张量σ的散度，并且
图 4.2.左：科学家兼工程师 C. Navier（1785 年 2 月 10 日至 1836 年 8 月 21 日）。右：G. Stokes 爵士（1819 年 8 月 13 日至 1903 年 2 月 1 日）。
118
数值建模简介
替换方程 （4.9） 中的结果，我们得到：
张量形式：ρa = F − ∇p + （λ+ μ） ∇（ ∇ · v ） + μv
成分形式：ρa = F−
∂p
∂x
+ (λ+ μ)
∂D
∂x
+ μ
∂v
∂ x
其中，拉普拉斯算子的向量为：v = ∇ ·∇v = ∇v
(4.10)
这种张量形式相当于一组分量形式的三个标量方程 （i = 1， 2， 3），称为纳维-斯托克斯方程。它们可以应用于描述任何介质中任何粘性流体的流动，包括储层岩石的孔隙网络。剪切动态粘度 μ 测量流体的流动阻力。在恒定温度下，液体的动态粘度会随着压力的增加而降低，但它对压力的依赖较小，而对流体温度的依赖更大。例如，对于 250C 和 15.0 MPa 的水，μ= 109 × 10Pa ·s;在相同温度和 4.0 MPa 下，μ= 106 × 10Pa ·s;而在 4.0 MPa 和 50C 的相同固定压力下，μ= 548 × 10Pa ·s （图 4.3）。只有当 p 达到较大值时，μon 压力的依赖性才会明显，例如在深度超过 2000 m 的储层中。
体动力粘度 λ 测量流体被压缩或膨胀的阻力。
实验发现，对于不同的流体 λ≥ 0 总是正确的 （Truesdell 1963）;λ= 0 仅当流体可压缩性完全可以忽略不计时，就像冷水一样。一般来说 λ≥ μ，商 λ/μ 在气体中较大，在稠密液体中较低。对于水，实验范围为 1 ≤ λ/μ≤ 3。对于气体，这个商要大得多;例如，大气条件下 CO 的 λ/μ≈ 10 （Truesdell 1963）。对于寒冷含水层中的水，在 25C 和 2.0 MPa 时，μ= 890 × 10Pa ·s 和 λ≈ 2600 × 10Pa ·s.对于地热系统，在 250C 和 20 MPa 时，液态水的剪切粘度为 μ= 110.4 × 10Pa ·s.在这些条件下，其体积动态粘度的相应值为 λ≈ 330 × 10Pa ·s.对于 300C 和 1.0 MPa 的蒸汽相，蒸汽的剪切粘度为 μ= 20.2 × 10Pa ·s;在相同温度和 8.5 MPa 下，μ= 19.7 × 10 Pa ·s. 图 4.3 显示了水的两相剪切粘度对温度的依赖性。请注意，在接近水的临界点（374.15C，22.12 MPa）时，两相的粘度趋于塌陷到共同值 47 × 10Pa ·s.Euler 方程适用于非粘性流体，因为它们忽略了流体粘度的影响。这些经典方程作为特例包含在 Navier-Stokes 系统中 （λ= μ= 0）。因此，Euler 方程的解只是实际流体流动问题的粗略近似。
0 40
饱和温度 °C）
动态粘度 （10
–6
帕 ·s)
液体粘度
蒸汽粘度
100
200
300
500
400
600
80
120
160
200
240
280
320
360
图 4.3.
剪切动粘度 μ 用于两相水。
268
数值建模简介
140
145
150
155
160
165
140
155
160
165
km
0
5
河 2
N
河 1
145
150
km
0
5
河 2
N
河 1
–1
–54
–15
–19
–17
–68
–55
–8
–6
–16
–35
–33
–19
–26
–7
–17
–41
–7
–2
–8
–39
4
19
19
42
58
71
19
44
83
47
29
3
1
27
17
23
15
30
22
11
26
45
24
30
1
3
9
30
44
35
12
254
243
244
134
240
141
189
98
172
132
145
130
105
45
11
45
73
65
48
100
48
121
190
167
69
20
36
37
41
39
41
22
44
51
53
66
65
63
68
45
39
14
18
49
63
49
5
1
50
6
5–6 6–7
7–8
3–4 4–5
>8
0–0.5 0.5–1
1–2
和观测值 （M）
2–3
模拟的差异
液压头：
5–6 6–7
7–8
3–4 4–5
>8
0–0.5 0.5–1 1–2
和观测值 （M）
2–3
模拟的差异
液压头：
山脉
山脉
模拟液压扬程 （m a.s.l.）
h
>
sim h
渗压计：h
–
sim h
h
<
sim h
(cm)
模拟液压扬程 （m a.s.l.）
h
>
sim h
渗压计：h
–
sim h
h
<
sim h
(cm)
d
c
图 6.9.
（续）
6.9.3
敏感性分析
敏感性分析是一种系统测试，用于评估模型参数的变化如何影响建模结果。模型的不同参数在特征范围内发生变化，该特征范围是根据参数在现实世界中行为的先验知识设置的。参数的可能变化由在现场观察到的参数的不确定性给出（另见第 6.12.1 节）。根据因变量（例如，水力水头、流体流速、溶质或热通量）观察到的模型响应的相对变化定义了影响最大的参数。模型更敏感的参数需要更多
270
数值建模简介
km
0
5
河 2
N
河 1
140
145
150
155
160
165
km
0
5
河 2
N
河 1
140
145
150
155
160
165
3–3.5 3.5–4
4–4.5
2–2.5
5–5.5
2.5–3
>5.5
4.5–5
0–0.5 0.5–1 1–1.5
由抽水引起 （m）
1.5–2
地下水位下降
3–3.5 3.5–4
4–4.5
2–2.5
5–5.5
2.5–3
>5.5
4.5–5
0–0.5
0.5–1 1–1.5
由抽水引起 （m）
1.5–2
地下水位下降
地下水流向井的路径
模拟液压扬程 （m a.s.l.）
地下水流向井的路径
模拟液压扬程 （m a.s.l.）
山脉
山脉
Well
垃圾填埋场
Well
垃圾填埋场
a
b
图 6.10.
饮用水井的捕获区，用于由最外层路径划定的井场的不同开采率;每 5 年一次的时间制定者：（a） Q= −0.4;（b） −1.0;（c） −0.7 米/秒。此外，井场开采导致的地下水位下降与自然初始地下水位的图。
干旱期将影响当地和区域水位。它可以预测污染物羽流在含水层中的传播，以及地热储层中的热传递。通过这种方式，优化地下水开采、修复受污染的含水层和优化地热田的开采将成为可能。
在该示例中，应调查是否可以在图 6.10 所示的地点建立井田。问题是来自垃圾填埋场的受污染水是否可以到达井场，以及为了
数值模型阐述过程
271
km
0
5
河 2
N
河 1
140
145
150
155
160
165
地下水流向井的路径
模拟液压扬程 （m a.s.l.）
3–3.5 3.5–4
4–4.5
2–2.5
5–5.5
2.5–3
>5.5
4.5–5
0–0.5 0.5–1 1–1.5
由抽水引起 （m）
1.5–2
地下水位下降
山脉
Well
垃圾填埋场
c
图 6.10.（续）
避免这种情况。此外，为了获得开采该油田的许可，因抽水而导致的天然地下水位下降不得超过井田周围半径为 3 公里的圆圈划定的区域外 3.5 m。
这项任务可以使用先前校准的数值模型来解决，并使用井田的不同开采率 Q 执行一系列数值模拟。图 6.10 显示了不同开采率的井田捕获区。对 φ= 0.08 的有效含水层孔隙度和 R = 1（保守示踪剂）的延迟因子进行了模拟。不应考虑污染物的横向分散。此外，地下水流路径上绘制的是指示旅行时间的时间标记（每 5 年）。根据 Q 变化生成的结果，估计最大值 Q = −0.7 m/s 是井场中允许的最大取水量，以便垃圾渗滤液不会进入其集水区。此外，第二个条件，即避免在 3 公里半径外出现超过 3.5 m 的洼地，也适用于该提款率。
6.11
模型有多好？评估不确定性
地下水系统的数值建模与许多不确定性有关。数值模型不确定性有两种类型，所选概念模型固有的不确定性和模型参数值的不确定性。数值模型所基于的概念模型是对真实场尺度情况的简化描述，其中包括许多简化的近似和假设，特别是关于含水层几何形状和模型参数的估计值、变量的敏感性、边界位置和初始条件及其空间和时间变化的近似和假设。
模型域边界是不确定性的来源：边界的位置可能会随时间变化，或者同一位置的边界条件（类型）可能会随时间变化。边界的条件（例如，特定的水力水头或溶质浓度、边界通量等）可能是时间的函数。边界磁通量难以测量，因此通常具有高度不确定性。参数值的空间分布可能在短距离内在多个数量级内发生变化。此行为是由异质性引起的
272
数值建模简介
在水力传导率 K、存储系数、S、有效孔隙度、φ 和分散性 α 上，这些值都是使用来自几个位置的离散值近似的，这些位置是在现场进行测量的。如果使用二维水平模型，则每个点（在水平平原中）的值还对应于含水层深度的平均值。其他参数值不仅取决于空间，还取决于时间（例如，水头、地下水补给、溶质浓度、温度），它们的不确定性是由于缺少时间数据测量，或者如果没有识别和考虑这些参数的时间依赖性，则仅使用时间平均值。不确定性还与流体和固体的汇和源有关，这可能是时间的函数，以及归因的传输机制和传输参数，包括描述化学反应、吸附过程和衰变过程的那些，这些还取决于固相组成、地球化学环境、温度等。
已经开发了不同的方法，将不确定性引入数值建模中，以解释输入参数的不确定性和概念模型的不确定性。对其预测结果中概念模型不确定性的评估通常是通过使用可用的不同模型进行模拟并评估为不同模型获得的预测范围来完成的（例如，Medina 和 Carrera 1996）。参数不确定性对建模结果影响的评估可以通过各种方法建立，这些方法可用于量化建模结果（预测）不确定性：线性近似、非线性近似和蒙特卡洛方法（参见 Carrera 等人，2005 年）。在后者中，各种可能性在大量模拟实现中表示，并获得结果参数分布的统计参数（例如，Carrera 等人，2005 年）。其他方法使用随机模型，其中各种系数表示为概率分布。
6.12
模型误用和错误
在数值建模的不同步骤中，可能会发生不同类型的错误和误用（参见 Mercer 和 Faust 1981，Mercer 1991）。它们可以分为四组：（1） 对要考虑的问题进行不正确的概念化，（2） 选择不合适的建模代码，（3） 不正确的模型应用，以及 （4） 对模型结果的误解。
如前所述，第一组与制定正确反映地下水流、溶质和热传输过程的准确概念模型有关，是制定合适的数值模型的先决条件。如果概念模型错误或不够准确，则获得的模拟结果不能反映自然系统的行为。在概念模型制定过程中，主要错误主要与模型域区域（含水层几何形状、边界位置、边界类型和边界值）的不当划定有关，对含水层的均匀性、各向同性和水力参数（以及它们各自的空间分布和随时间变化）的错误假设有关，对发生的运输过程的错误假设， 以及维度选择不合适（例如，在需要 3D 模型的情况下使用 2D 模型）。
对于第二组，通常可以观察到建模者在不需要此类模型的情况下使用高度复杂的模型程序，因为没有足够的现场数据来支持它或因为目标不需要它。
如果使用不正确的输入数据，并且没有正确选择网格或网格大小以及时间步长的间隔，则会出现不正确的模型应用程序。此外，有时所选的模型代码与所选的概念模型不兼容。在模型校准期间，选择不合适的校准参数、校准周期以及使用在与建模时间间隔不同的条件下校准的模型都会产生额外的错误。这些错误可能导致错误的结果和对建模结果的错误解释。
数值模型阐述过程
273
6.13
模型构建示例 - 含水层污染评估
在下文中，我们提供了一个理想化的示例，其中我们展示了建模方法的逐步过程。必须考虑到，为了不分散对主要建模目标的注意力，我们进行了几次简化。
6.13.1
情况和任务
在拉斯帕尔马斯市附近（15,000 名居民），地下水是从位于不同物业的井中提取的。这种水主要用作饮用水，其次用于灌溉（图 6.11）。两年前，一家于 1989 年开始生产的塑料材料工厂建在该镇的西部边界（作为时间参考，“今天”对应于 1991 年 6 月）。该镇本身的供水来自位于东南 3 公里处的自来水厂（图 6.11）。该自来水厂以 Q = −0.05 m /s 的恒定速率从井中抽水。根据一些农民的投诉，他们注意到水质与往年相比有所下降，1991 年 6 月 16 日，镇政府委托在几口井中进行了取样和化学分析活动。
结果显示，与邻近地区的钠钠（最大 750 mg/l）、氯化物 Cl（最大 1030 mg/l）、硫酸盐 SO 4 （1348 mg/l） 和羟基硼酸根离子 B（OH）（最大 17.5 mg B/l） 的浓度相比，值显著增加（表 6.6，图 6.11a）。水质的下降一方面是由于它的咸味，另一方面是由于高矿化和羟基硼酸根离子的浓度，这阻止了地下水用于人类消费和灌溉，就像迄今为止使用的那样，或者允许它的使用程度要小得多。
W22
W24
W16
W15
W14
W19
W25
W31
W34
W30
W33
W26
W36
W32
W35
W29
W28
W23
W20
W27
W21
W17
W18
W12
W11
W13
W7
W9
W6
W10
W5
W8
W1
W4
W2
P5
PW
P8
P7
P6
P4
P3
P2
P1
W3
152
153
500
0
m
a
Las
帕尔马斯
P
166
164
162
160
158
156
154
167
165
163
161
159
157
155
151
Las
帕尔马斯
P
b
500
0
m
> 100 毫克/升 Cl
> 1 毫克 B/l
测压高度（观测）（m a.s.l.）
–
图 6.11.拉斯帕尔马斯的现场示例：（a） 钻探井 W1 至 的研究区域位置
W26，渗压计，P1 到 P7，抽水井，PW，和塑料材料工业厂房，P。此外，还显示了氯化物含量超过 100 mg/l 和硼含量超过 1 mg/l 的区域（以地下水中发现的羟基硼酸根阴离子的形式存在）;（b） 等高线图显示了 1990 年 6 月至 1991 年 6 月期间观测到的地下水位的平均值，在自来水厂的 W1 至 W36 井、渗压计 P1 至 P7 和抽水井 PW 处测量（数据见表 6.6）。
274
数值建模简介
表 6.6.
拉斯帕尔马斯的现场示例：地下水中的氯化物、硫酸根和羟基硼酸盐离子浓度。样本于 1991 年 6 月 16 日采集，平均渗压测量值是 1990 年 6 月至 1991 年 6 月期间的测量值（1990 年 6 月 11 日、1990 年 9 月 15 日、1990 年 12 月 1 日和 1991 年 6 月 16 日的测量）。
采样
测压
Cl
SO
B（羟基）
点
级别 （M A.S.L.）
（毫克/升）
（毫克/升）
（毫克 B/l）
W1
164.6
1.3
21.5
<0.1
W2
166.8
1.8
25.0
<0.1
W3
167.6
0.7
19.7
<0.1
W4
165.7
1.4
13.7
<0.1
W5
164.0
1.5
17.3
<0.1
W6
161.8
1.1
22.5
<0.1
W7
166.0
0.9
21.7
<0.1
W8
162.7
1.6
19.6
<0.1
W9
159.2
1.8
24.8
<0.1
W10
161.9
1.2
21.9
<0.1
W11
164.3
0.8
19.1
<0.1
W12
162.2
2.8
56.5
<0.1
W13
163.4
1.0
25.2
<0.1
W14
159.0
1.7
23.6
<0.1
W15
157.7
0.5
21.7
<0.1
W16
156.1
1.3
20.5
<0.1
W17
161.8
1.5
21.9
<0.1
W18
160.4
105.0
267
0.1
W19
157.6
1.5
23.4
<0.1
W20
162.9
1.4
26.0
<0.1
W21
158.6
1030.0
1348.0
<0.1
W22
162.0
780.0
1037.0
17.5
W23
162.0
0.9
23.7
<0.1
W24
153.5
1.8
19.5
<0.1
W25
154.4
1.3
17.8
<0.1
W26
156.8
1.3
14.9
<0.1
W27
158.8
1.8
21.1
<0.1
W28
160.5
1.7
23.7
<0.1
W29
159.2
1.6
25.0
<0.1
W30
153.7
1.1
24.8
<0.1
W31
132.6
1.5
12.9
<0.1
W32
156.6
1.3
21.8
<0.1
W33
154.5
1.3
25.7
<0.1
W34
151.8
1.6
21.6
<0.1
W35
157.4
0.9
15.8
<0.1
W36
154.2
1.6
23.7
<0.1
P1
−
1.0
19.8
<0.1
P2
−
1.5
17.7
<0.1
P3
−
10.5
56.8
0.4
P4
−
1.4
23.8
<0.1
P5
−
1.6
23.9
<0.1
P6
−
0.9
21.3
<0.1
P7
−
1.7
19.7
<0.1
PW
150.9
194
245.7
<0.1
在这种情况下，存在许多问题：
•
谁造成污染？最近建造的工业设施（可能是一个污染源）是否要对这种情况负责？这种责任可以证明吗？如果工厂不负责，污染物是如何进入系统的？
数值模型阐述过程
275
•
自来水厂是否存在严重或潜在风险？如果有，未来污染将如何演变？污染羽流何时到达抽水井，何时（如果有的话）会消失？这些问题必须从两个角度来回答：一方面，考虑最佳情况，即立即关闭污染源;另一方面，最坏的情况是污染过程继续保持不变。对于这两种情况，都应该确定将到达抽水井的浓度，包括其最大值的大小以及它们在时间上的分布。必须研究明显污染的区域，以及预计会发生污染的区域。农场水井必须执行相同的程序。
•
如有必要，是否有可能或合理地采取保护或卫生措施来维持现有自来水厂对拉斯帕尔马斯的供水？哪些解决方案可用于向农场供水？
6.13.1.1 现有字段数据和信息
关于非承压含水层的范围和岩性特征的唯一可用信息来自 7 个渗压计（P1 至 P7）和一口抽水井 （PW），其中有地质钻探剖面。它们都到达了第四纪含水层的底部，该含水层由不透水的块状石灰岩组成。含水层由中等粒粒的沙子和砾石组成，淤泥百分比低于 5%。在 W1 至 W36 的 36 口钻井中，可以测量渗压液位，并在 1990 年 6 月至 1991 年 6 月期间进行了为期一年的系统测量。平均值见表 6.6。1991 年 6 月 16 日，采集样本并进行相应的分析。表 6.6 列出了 Cl、SO 和 B（OH）（如 B）获得的浓度值。
6.13.2
调查计划的设计
如果模型用于处理污染事故，建议遵循一定的顺序，例如 DVWK （1989） 编制的顺序：
•
初步研究 — 污染风险评估
•
评估方法 - 水力水头测量和污染物浓度测定
•
简单的模型 - 理想的系统
•
流动时间计算 - 纯平流输运建模作为第一种方法
•
考虑分散、吸附和其他过程的输运建模
6.13.3 初步调查 - 污染评估
在拉斯帕尔马斯和自来水厂之间观察到的氯化物浓度高达 1030 毫克/升，而在抽水井处观察到的氯化物浓度达到 194 毫克/升（阈值为 250 毫克/升），因此污染风险可以归类为急性。因此，必须进行研究活动，以确定污染物的真实分布及其未来的传播。
6.13.4
获取地下水位和污染物浓度
可以通过以下方式对污染物传播进行初步粗略评估：
•
评估地下水流的速度和方向
•
使用已经存在的污染物
在第一种情况下，进行地下水位测量或气压测量，而在第二种情况下，使用化学分析中获得的物质浓度值。对于我们的示例，可以使用现有的渗压测量数据，以及要研究的每个离子的分析结果。根据测压数据，可以绘制等高线图（图 6.11b）;
276
数值建模简介
这张地图显示了从大约 WNW 到 ESE 的地下水流。通过观察溶质浓度的分布，可以表明对应于不同物质的污染羽流具有相似的形状，其来源位于工业厂房（图 6.11a）。另一方面，在可能成为另一个潜在污染源的城镇与可疑的工业厂房之间，与自然背景值相对应的浓度较低。因此，工业厂房似乎是唯一的污染源。因此，污染源是局部的，其开始时间是已知的，即 1989 年 6 月。然后，自来水厂直接位于点源的下游。其迁移率（延迟因子 R ≈ 1）大于其他溶质的氯化物污染羽流已经达到它（图 6.11a）。其 Clconcentration 为 194 mg/l，仍低于 250 mg/l 的阈值和约 350 mg/l 的品尝阈值。污染源点与自来水厂之间的地下水流平均线速度（平均孔隙速度）v 可以根据水力梯度 i≈ 0.0072、含水层介质（含细砾石的中沙）的近似水力传导率 K（1 × 10≤ K ≤ 1 × 10m/s）和有效孔隙度 φ（0.05 ≤ φ≤ 0.2）进行近似估计。使用 v= Ki /φ，v × 3.6 10 和 1.4 × 10m/s 之间的值。
从这些信息可以得出以下结论：
•
可以确定导致地下水污染的责任方，从而确定污染点的来源。证实了污染是由最近定居的塑料材料工业工厂引起的怀疑。
•
从源头开始，污染羽流沿 WNW-ESE 方向前进，氯化物的传播速度可以归类为地下水的传播速度，为 3.6 × 10≤ v≤ 1.4 × 10m/s。
•
受影响或可能受影响的区域是划定的。
6.13.5
地下水流和平流输运作为第一种方法
到目前为止进行的研究活动只是一般的近似值;为了更精确地描述含水层中的运输过程，需要研究和定量描述水力状况。这对于描述平流传输至关重要，它经常通过使用解析模型（在我们的例子中是数值模型）在含水层中的物质传播中起着主导作用。地下水流速模式的确定以及其他水文地质参数特别值得关注。使用简单的数值模型，例如 ASM 建模代码（参见第 6.8.3 节），可以模拟地下水流，并可以根据污染源点和从该点开始的地下水流路径计算水颗粒的流动时间。
6.13.5.1 划分建模区域和含水层几何形状
为了简单描述模型边界，将要建模的污染区放置在其两个边界与地下水位等值线图（西北和东南边界）平行的位置，而另外两个边界大致垂直方向，这意味着没有地下水流通过这些边界（不透水的 NE 和 SW 边界）（图 6.12a）。然后，水力水头在 NW 和 SE 边界的边界条件对应于相邻等值线的边界条件（图 6.12a）。根据在钻孔处采集的数据，绘制了含水层底部的等深图，用作模型的垂直下限（图 6.12b）。
6.13.5.2 水文地质参数
1990 年 6 月至 1991 年 6 月期间在不同季节进行的水力水头（测压水位）测量（1990 年 6 月 11 日、1990 年 9 月 15 日、1990 年 12 月 1 日和 1991 年 6 月 16 日的测量）表明，水位的季节性变化低于 0.5 m，而且没有
数值模型阐述过程
277
Las
帕尔马斯
P
153
P
帕尔马斯
Las
a
b
c
d
500
0
m
500
0
m
152
153
166
164
162
160
158
156
154
167
165
163
161
159
157
155
151
152
153
166
164
162
160
158
156
154
165
163
161
159
157
155
151
建 模
每 100 天带标记的流径
观察
压压高度：
148
146
144
Las
帕尔马斯
P
142
144
146
142
144
146
148
500
0
m
含水层基础 （m a.s.l.）
Las
帕尔马斯
P
500
0
m
152
153
166
164
162
160
158
156
154
167
165
163
161
159
157
155
151
模型区域
恒定液压头单元
W
图 6.12.拉斯帕尔马斯的现场示例：（a） 含水层底部的位置，由
来自压力计 P1 至 P7 和 PW 孔的基底深度数据;（b） 研究和建模区域的位置，具有有限差分模型单元;（c） 在建模区域中以数字方式建模的地下水位等高线图;（d） 模拟建模区域中的地下水流路径，每 100 天按时间顺序划分一次，井的集水区（深灰色）。图 （b） 到 （d） 还显示了基于现场数据位于建模区域之外的地下水等值线图。
对地下水流速模式有重大影响。因此，数值模拟是在假设地下水流静止的情况下，使用测压水平的平均值来完成的（表 6.7）。为了确定水力特性，在压力计 P1 到 P7 上进行了泵送测试，从中获得了水力传导率 K 和存储系数 S（表 6.8）。K 值范围从 3.5 × 10 到 6 × 10m/s，在目前的含水层中对应于有效孔隙度的 S 值在 0.08 到 0.14 之间变化，并且没有显示出明显的空间模式。但是，请记住，这些值基于离散点的少量现场测量。在实践中，需要在建模区域内外获得更多的测量值。这将允许建模者划定值范围，其中可以预期离散值。这与灵敏度分析一起，用于评估特定参数值的变化对仿真结果的影响程度。
278
数值建模简介
表 6.7.拉斯帕尔马斯的现场示例：用于压力计 P1 至 P7 和抽水井 PW 的含水层底部的海平面以上高度。
炮眼/
渗压计
P1
P2
P3
P4
P5
P6
P7
PW
含水层基地
147.7
143.7
142.3
145.0
148.1
145.1
143.2
144.2
（M A.S.L.）
表 6.8.
拉斯帕尔马斯的现场示例：水力传导率 K 和存储系数 S，通过压力计 P1 至 P7 的泵送测试获得;含水层厚度 b ，从测压水平与含水层底部之间的差异获得。
嗯，压力计
P1
P2
P3
P4
P5
P6
P7
含水层底部（地表以下 m）
147.7
143.7
142.3
145.0
148.1
145.1
143.2
测压级 h （m a.s.l.）
161.8
165.0
164.0
164.4
155.4
160.1
156.1
含水层厚度 b（m）
14.1
21.3
21.7
19.4
7.3
15.0
12.9
水力。电导率 K（10m/s）
5.2
6.0
5.0
3.5
4.5
6.0
4.9
存储系数 S （ad）
0.10
0.14
0.11
0.08
0.10
0.11
0.08
可用于获得一系列可能的解。降雨对地下水的补给被忽视了。
6.13.5.3 地下水流线和流速模拟
使用上述数据和信息，可以构建数值流动模型。校准后，使用 K 作为校准参数，将模拟的压能级与现场观察到的能级相匹配，即可进行模拟。
6.13.5.4 结果 地下水流速模拟获得的结果仅限于溶质迁移和污染物迁移的解释，但给出了重要的初步结果。原因是没有考虑溶质分散效应和吸附、解吸和分解过程。然而，在氯化物作为污染物的情况下，可以假设它不受吸附过程的影响，并且表现为理想的示踪剂。相比之下，羟基硼酸根离子经历吸附过程，这就解释了为什么氯化物羽流已经到达饮用水井 PW，而硼浓度仍然对应于背景值（表 6.6）。
基于我们的无色散纯平流污染物传输方法作为我们问题的第一个近似值，我们从数值建模中获得了以下结果：
•
根据流线和流线上的时间标记，从污染点源（100 m 宽）离开的水颗粒需要 500 到 600 天才能到达自来水厂的抽水井。表现为理想示踪剂的物质（如 Cl）也会发生同样的情况。
•
抽水井捕获区的宽度可以从图 6.12d 中确定为 W = 700 m。这意味着该井接收来自 700 m 宽波段的水，其中包括污染源和整个羽流（图 6.12d）。根据氯化物流入量 C= 1000 kg/天 （11,574 mg/s） 和羟基硼酸盐阴离子流入量 C= 50 kg B/天 （579 mg B/s） 的估计值，并使用生产井的抽水速率 Q= −0.05 m/s （50 m/s），可以估计抽水井的最大预期浓度。从该计算中，获得了氯化物 C= C/Q = 231 mg/l 和羟基硼酸根阴离子 C= C /Q = 12 mg B/l 的最大预期浓度。
数值模型阐述过程
279
•
最高氯化物浓度 关注 略低于 WHO 和欧盟承认的阈值（分别为 200 和 250 毫克/升），而最大硼浓度 关注 远高于 EPA 设定的阈值，不超过 1 毫克/升硼（WHO：考虑的一年为 0.1 毫克 B/l）。
•
受当前和未来污染影响的区域已被划定;它从污染点源扩散到饮用水井，在那里通过抽水井从含水层中提取许多污染颗粒，从而防止污染超出该位置。
•
如果未消除污染点源，则氯离子和羟基硼酸盐阴离子浓度将出现准平稳分布。
•
如果污染立即中断，则在氯化物（理想示踪剂）的情况下，污染将在 t = 500 至 600 天内结束。
6.13.6
含分散、吸附和结果解的输运模型
6.13.6.1 简介
考虑分散、吸附和其他过程的传输模型的制定需要更多的信息（参见第 6.5 节和第 6.7 节关于数据需求和数据收集）以及更密集的采样网和更频繁的采样。因此，与以前的研究相比，它在时间和金钱方面要困难得多，成本也更高。必须仔细评估每个特定情况，以确定从第 6.13.5 节中构建的纯平流模型获得的结果是否已经足以给出问题的所有答案。如果需要考虑吸附，我们需要知道是否可以从文献中获取溶质特异性吸附数据， 或者是否必须在实验室进行吸附实验。这同样适用于获取纵向和水平色散率值。
在拉斯帕尔马斯的例子中，地下水流场在之前的研究中已经非常有名（见 6.13.5），因此满足了制定溶质输运模型的基本条件之一。相反，关于溶质浓度的信息非常有限（对于构建传输模型来说太有限了），因为采样和分析只能进行一次，而且污染羽流覆盖的区域只有 6 个采样点。没有与含水层的分散参数对应的数据。关于污染物与含水层材料以及溶解在地下水中的其他物质之间的交换过程，只有参考书目可用。据此，氯化物可用作理想的示踪剂，而硼在目前的 pH 值下以羟基硼酸盐阴离子的形式存在，经历吸附和解吸过程，主要吸附到粘土矿物、铁、氢氧化锰和氧化物以及有机物质上。
基于这种有限的数据情景，我们研究的继续有两种基本可能性：
•
如果可能，可以根据第 6.5 节和第 6.6 节中关于数据需求和数据收集的内容进行详细调查，以便在给定的时间范围内获取缺失的数据。通常，除特殊情况外，传播模型只能在经过长时间广泛研究的领域进行详细说明，因此有丰富的数据可用于校准模型。
•
因此，在实践中，通常会采用折衷解决方案，并且通过估计或假设（例如，关于吸附参数和含水层分散性）来弥补数据的缺乏。在这些情况下，可以对要描述的含水层进行敏感性分析，改变不同的参数，以确定这些参数对系统行为的影响。自然，这种方法降低了所获结果的可信度。为了产生可靠的定量预测，应该能够使用模型验证现场测量，包括流型、浓度分布及其随时间的变化。
280
数值建模简介
在本例中，使用了 6.13.5 节中的稳态地下水流模型。此外，还执行了两次示踪剂测试，以获得纵向和横向色散率的值（有关示踪剂测试，请参见第 6.7.2.5 节）。同样，我们需要记住，这些含水层参数是根据仅两个离散点的现场测量确定的，因此，即使结果是正确的，它们也可能仅适用于被测试的小局部区域。在实践中，需要更多的确定来确定值的范围，这将使建模者能够获得一系列建模结果。对于每项示踪剂测试，使用两个孔/渗压计，一个作为进样孔，另一个作为采样点。如果不钻探新的渗压计，测试仅限于几个潜在位置，因为大多数井之间的距离大于地下水颗粒传播时间的 100 天。其中一项测试在 W5 和 W10 井之间进行，而第二项测试在 W17 和 W28 井之间进行（图 6.11）。在这两种情况下，相应的第一口井作为注入井，第二口井作为采样点。这些示踪剂测试的评估分别产生有效孔隙率 φ 以及纵向和横向分散率 α 和 α 的值（表 6.9）。通过示踪测试获得的数据对建模区域中的两个离散点有效，在解释建模结果时必须考虑到这一点。 由于没有关于含水层基质中吸附材料和地球化学条件的信息，并且由于时间和资源不允许收集这些信息，因此羟基硼酸盐阴离子的吸附行为应作为一个非常粗略的近似值进行评估。这个近似值是通过比较硼羽流和氯化物羽流来获得的，氯化物羽流可以被认为是理想的示踪剂，并且不会发生吸附和延迟。利用这些有限的数据，我们对两种可能的情况进行污染物迁移模拟，分别对应于最佳和最坏情况：
•
污染源使含水层的污染保持不变。
•
污染源立即关闭。
6.13.6.2 模拟污染物永久流入的溶质传播
第 6.13.5 节中阐述的校准数值地下水流模型现在正在扩展到考虑扩散和吸附过程的传输模型。但是，此模型的准确性有限，在解释仿真结果时需要考虑这一点。将建模的溶质分布与观察到的溶质分布相匹配的校准并不准确，需要进行一系列浓度测量才能达到效果。使用从示踪剂测试中获得的色散率的平均值，对 Clpropagation 的模拟产生了时间 t = 730 天的污染羽流的 Clconcentration 分布，该分布与在现场观察到的分布相匹配（图 6.13a）。另一方面，与 Clpropagation 相比，硼的传播延迟。这表明羟基硼酸盐阴离子正在被含水层介质吸附。使用羟基硼酸盐阴离子的延迟因子 R 校准模型，并假设线性吸附，当值 R = 2.5 时，模拟的污染羽流与观察到的污染羽流最匹配（图 6.13b）。
考虑到上述假设，模拟了一些受污染井的浓度随时间的变化。氯化物浓度曲线显示，位于污染羽流边缘的受污染井 W22、W21、W18、PW 和 P3
表 6.9.拉斯帕尔马斯的现场示例：通过示踪剂测试获得的有效孔隙率和分散性值。
示踪剂 I 的进样井
W5
W17
采样观察井 U
W10
W28
有效孔隙率 φ（ad）
0.12
0.095
纵向分散率 α（m）
12
9
横向分散率 α（m）
1.0
0.7
数值模型阐述过程
281
W28
W28
W17
1
10
30
100
500
1000
W22
W15
W14
W19
W25
W33
W26
W27
W21
W17
W18
W12
W11
W13
W10
W5
W8
W4
P5
PW
P7
P4
P3
P2
Las
帕尔马斯
P
P
帕尔马斯
Las
P2
P3
P4
P7
PW
P5
W4
W8
W5
W10
W13
W11
W12
W18
W21
W27
W26 W33
W25
W19
W14
W15
W22
a
b
500
0
m
500
0
m
Cl
B（羟基）
–
720 天和 2000 天后的分布 （mg/L）
（mg B/l） 之后
720 和
2000 天
0 图 6.13.来自拉斯帕尔马斯的现场示例：1991 年 6 月 16 日（=730 天）和污染物输入开始后 2000 天的污染物分布建模：（a） 氯离子，（b） 羟基硼酸根离子。730 天和 2000 天的 Cldistributions 相同。
days
年
0
1000
800
600
400
200
0
40
30
20
10
PW
W21
W22
W18
P3
PW
P3
W18
W21
W22
0
0
400
800
1200
1600
3
2
5
4
1
氯化钙 （mg/l）
B（OH） （毫克 B/l）
days
年
0
0
400
800
1200
1600
3
2
5
4
1
图 6.14.溶质浓度的演变，针对受影响的扇区建模：（a） 氯化物，（b） 羟基 -
硼酸盐阴离子（值为 B）（t = 0 年对应于 1989 年 6 月污染物流入的开始;t = 2 年对应于现在的时间，即 1991 年 6 月）。
282
数值建模简介
目前（t = 2 年）处于稳定阶段，这意味着氯化物浓度在未来不会再变化（图 6.14）。另一方面，对应于硼的浓度曲线表明，这些井尚未达到稳态，浓度在未来会增加（图 6.14）。
因此，可以得出以下结论：
•
只有 W22、W21、W18 和 PW 孔的氯化物浓度为 860、980、135 和 180 mg/l，将来不会改变。W22 和 W21 井无法使用，其余井（包括抽水井）低于阈值。
•
位于受影响区域内的 W22、W21、W18 和 PW 井的羟基硼酸根离子浓度分别达到 20、0、0.5 和 0.00 mg /l 的值（时间 t = 2 年）。这些值将在未来增加，直到达到 49、41、10 和 9.3 mg B/l，并且许多井将变得毫无用处。在水厂的抽水井将在 t = 3.5 年时超过 1 mg B/l 的阈值，即从现在起 1.5 年。W22 井由于氯化物浓度高而无法使用，已经超过了硼的最大值。井 W21 和 W18 将在大约 7 个月内达到此阈值，即 t = 2.6 年。
6.13.6.3
模拟溶质传播，立即暂停污染物流入
如果立即关闭污染源（t = 730 天），则获得以下结果：
•
W22、W21、W18 和 PW 孔的氯化物浓度为 860、980、135 和 180 mg/l，在 0.5、0.9、1.0 和 1.5 年内大致保持不变。在 1.1、1.4、1.2 和 1.7 年后，将达到 Clbackground 值。对于 W22 和 W21 孔，氯化物浓度分别需要 0.8 年和 1.2 年才能低于可接受的阈值水平;在此之后，这两个井将再次可用。
•
由于羟基硼酸根离子的延迟作用，高浓度井 W22、W21、W18 和 PW 分别为 2.8、4.0、3.7 和 4.5 mg B/l，需要比氯化物更长的时间才能低于 1 mg B/l 的允许阈值。
在评估这些结果、污染物浓度及其时间变化时，必须考虑到无法对传输模型进行精确校准，因为浓度测量值仅存在于一个离散时间。因此，将来需要收集新的浓度数据，以便与通过建模获得的结果进行比较。在出现差异的情况下，模型必须适应新的测量数据，并且必须相应地调整结果和预测。
6.13.6.4 自来水工程：诊断和推荐的解决方案
抽水井仅受到硼污染的威胁，在浓度超过阈值水平之前可以再使用 1.5 年。如果工业厂房立即停止污染，该井将在 9 mg/l 的浓度下再使用 3 年。可以采取不同的措施来提交给统治组织。在成本评估研究之后，他们将选择最合适的选项或多个选项的组合：
•
必须立即关闭污染源。建议将有害物质浓度高的生产液体流出物（10 米/天）转移到底部密封的蒸发池中。液体蒸发后，应将剩余的固体存放在合适的储存库中。在建造这些水池时，污水可以储存在水箱或合适的沉积物中。
•
自来水厂的水必须经过处理以降低硼浓度，预计硼浓度将超过最大阈值 3 年，以使水可饮用。
数值模型阐述过程
283
•
此外，可以在污染羽流内钻一口卫生井，以提取和处理受污染的水。
•
自来水厂可以在污染羽流之外但在目前使用的抽水井附近钻一口新的抽水井。如果使用此选项，则新位置将由详细的仿真模型确定。这种可能性具有以下优点：（1） 费用应由污染行业支付;（2） 4.5 年后（当旧井再次适合时），拉斯帕尔马斯镇将增加 2 口用于饮用水供应的水井，鉴于居民不断增长的需求，这一点很重要。然而，需要考虑到的是，如果目前正在使用的井关闭，那么污染羽流将进一步向下游延伸。因此，必须进行研究以确定位于该井下游的其他井是否存在污染风险。如果是这样，建议井继续运行，以便通过去除所有污染物来清洁受污染的含水层。
6.13.6.5 农场水井：诊断和推荐的解决方案
所有农场井中只有 W1、W22、W11 和 W18 井不可用或将不可用。由于这些农场需要的水很少，因此可以在污染持续期间通过油罐车供水。灌溉所需的水可以由位于邻近农场的水井提供（有相应的货币补偿），并通过临时灌溉渠道运输。第二种选择是暂停灌溉并向对污染负责的一方提出索赔，以弥补农业产量减少或不足。
第 7 章
参数识别和逆问题
Longina Castellanos 和 Angel Pérez
“逆问题在地球物理学中起着至关重要的作用，因为该领域的主要任务之一是探测地球内部，这既是出于经济原因，如石油勘探，也是为了追求有关我们星球的学术知识。”
Roel Sniedery 等人 （1998）
7.1
介绍
研究含水层系统或任何一般物理系统的科学程序可以分为三个部分：
•
系统参数化：发现一组最小的模型参数，其值完全表征系统。
•
直接问题或直接建模：发现物理定律（例如，地下水流方程，见第 4 章），允许对给定的参数值（例如，含水层的水力传导率，参见第 4.3 节）进行预测。
•
逆向建模：使用观察到的参数的实际测量值来推断模型参数的值（例如，水力传导率）。这个推理问题被称为逆问题。
在第 5 章中，我们研究了获得直接问题的数值解的方法。例如，这些是有限元法、线法和有限差分法。在本章中，我们从微分方程的解测量中研究参数识别或重建未知系数的逆问题。
在实践中，通常可以使用代表所讨论现象不同状态的一般行为的概念模型（参见第 6.3 章和第 6.4 章）。让我们考虑一下简单的微分模型：
dφ
dt
= −xe+ xCos（t）
(7.1)
假设参数 x 和 x的值是已知的。然后，可以使用一些适当的数值方法获得方程 （7.1） 的近似解，从而解决直接问题 （DP）：我们知道原因（参数 x 和 x）并评估影响（φ（t） 的数值），这使我们能够了解状态行为随时间的变化。
假设我们不知道参数的值，但有可能获得 φ（t） 值在 m 个不同时间的实验测量值（数据）。目标是使用这些数据来估计最合适的模型中的参数 x 和 x 的值。这是一个逆问题 （IP）：我们知道效果 （φ（t）， i = 1， . . . ， m） 并寻找原因 （xand x）。
∗ 资料来源：Roel Sniedery、Malcolm Sambridgez 和 Fernando Sansó：地球物理学中的逆问题：缩小理论与实践之间的差距。客座编辑序言，Inverse Problems 14 （1998）。
285
286
数值建模简介
P = 参数空间
–1
–1
O = 观测空间 = 直接问题 = 逆问题
影响
原因
O
P
图 7.1.
反直接问题的图形表示：原因包括：例如，水力传导率;效果是，例如液压头。
我们首先简单地定义逆问题的解决方案，即使用效应的测量值来确定原因。这种观点与相应的直接问题形成鲜明对比，后者的解决方案涉及根据其原因的完整描述来寻找效果。图 7.1 显示了逆问题和直接问题之间关系的图形方法。
这些问题通常分为以下几个术语：
•
识别或重建，寻找观察到的效果的原因。
•
控制或设计，寻找预期效果的“一个”可能的原因。
这两个问题是相关的，但由于目标不同，也存在一些数学后果。在识别问题中，理想的属性是解决方案的唯一性（可识别性），理想情况下，观察到的效果有一个特定的原因，人们希望获得该结果。在控制或设计问题中，唯一性并不重要，因为非唯一性仅意味着设计目标可以通过不同的策略来实现，因此，人们有额外的自由度（例如，纳入进一步的设计目标）。
逆问题可分为：
•
连续：大多数逆问题都属于此类型。要估计的参数集是多个变量中的连续函数。例如，含水层中的导水率是空间变量的连续函数。
•
Discrete：需要估计的模型参数数量有限，例如，上面示例中的参数 x和 x。
有时，问题本身在本质上（现实世界）是连续的，但由于计算原因而离散化。例如，使用在含水层区域设计的网格对含水层中用于查找水力水头的偏微分方程 （PDE） 进行离散化，从而产生离散方程组。
逆问题的模型函数，可能是微分方程组的解，如例 （7.1） 所示，可能线性或非线性地取决于参数。在方程 （7.1） 中，微分模型的解析解为：
φ（x， x; t） = xe+ xSin（t）
(7.2)
假设我们有 m 个数据点 yi 表示 φ 在不同时间的行为，我们需要找到 x和 x，使得每个 yi 到相应的模型函数 φ（t） 的距离对于每个 t，i = 1， . . . ， m，最小（见图 7.2）。在这个例子中，直接问题 [模型函数 （7.2）] 的解线性取决于参数 x 和 x，因此我们需要求解线性逆问题。
参数识别和逆问题
287
t
( )
t
t
t
t
t
t
y
图 7.2.模型 φ（t） 和观测数据 yi .
另一方面，如果差分模型为：
dφ
dt
= −xe+ xCos（xt）
(7.3)
则 model 函数为：
φ（x， x; t） = e+ Sin（xt）
(7.4)
因此，直接问题 [模型函数 （7.4）] 的解非线性地取决于参数 x 和 x，我们需要求解一个非线性逆问题。线性和非线性情况的目标都是使每个 t 处的计算解，即 φ（x， x; t ） = y 尽可能接近观测数据 yi ，对于所有 i = 1， . . . ， m。
在实践中，最常见的情况是微分模型的解析解不可用，而 DP 是通过数值求解的，因此很难确定 IP 是线性的还是非线性的。在微分方程上线性出现的参数并不意味着它在解中也是线性的。例如，方程式：
dϕ
dt
= −xφ
是 x 的线性函数。然而，它的解 φ（t） = e） 是参数 x 的非线性函数。
测量最近解的一种方法是定义某种距离。
最常见的原因之一是数据点和计算值之间的平方距离，原因有几个。在这种特殊情况下，希望获得值 x 和 x 对，使得总和：
m
∑
i=1
（φ（x， x; t ） − yi ）=
m
∑
i=1
(y
− 易 ）
(7.5)
最小化，并且参数 x 和 x将是达到总和最小值的参数。
一般的数学公式是所谓的最小二乘问题 （LSP），它如下（参数在模型函数中出现的特定形式使