地下水和地热系统数值建模简介
质量、能量和溶质的基本原理
多孔弹性岩石中的运输
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
•
自来水厂是否存在严重或潜在风险？如果有，未来污染将如何演变？污染羽流何时到达抽水井，何时（如果有的话）会消失？这些问题必须从两个角度来回答：一方面，考虑最佳情况，即立即关闭污染源;另一方面，最坏的情况是污染过程继续保持不变。对于这两种情况，都应该确定将到达抽水井的浓度，包括其最大值的大小以及它们在时间上的分布。必须研究明显污染的区域，以及预计会发生污染的区域。农场水井必须执行相同的程序。
•