Page 114 - 2025年第56卷第1期

基于中国地质调查局地质云数据库的地质钻孔数据， 并结合历史勘探资料， 深入分析了研究区域
的水文地质条件， 进而利用地下水模型系统（Ｇｒｏｕｎｄｗａｔｅｒ Ｍｏｄｅｌｉｎｇ Ｓｙｓｔｅｍ， ＧＭＳ）将含水层概化为三层
非均质潜水含水层， 水流概化为考虑变密度因素的三维非稳定流。
２．３．２ 边界条件概化及源汇项 研究区北部边界为莱州湾， 概化为海边已知水头边界； 库区东边为内
陆低山丘陵区， 由下第三系泥岩、 砂质泥岩组成； 库区西边为剥蚀台地， 下部为下第三系的泥岩、 泥
灰岩， 可概化为定流量边界， 但单位吸水量相对较小。 库区整体由南向北倾斜， 南部山丘区分布有下
元古界的片岩等， 补给性较好， 概化为定流量边界。 含水层顶部边界为潜水面， 底部边界概化为含水
层隔水底板。 在溶质运移模型中， 研究区北部为海洋边界， 概化为定浓度边界， 其余边界概化为已知
对流－弥散通量边界； 含水层顶部概化为零通量的水动力弥散通量边界， 底部隔水边界与外界没有溶
质的交换。 此外， 研究区内有黄水河经过， 主要的源汇项包括大气降水入渗、 河道补给、 定水头边界
补给、 河道排泄、 地下水开采、 边界侧向流出等。 研究区的黄水河采用 ＭＯＤＦＬＯＷ 中的通用水头边界
（Ｇｅｎｅｒａｌ Ｈｅａｄ Ｂｏｕｎｄａｒｙ， ＧＨＢ）子程序包处理， 地下坝采用水平流动障碍（Ｈｏｒｉｚｏｎｔａｌ Ｆｌｏｗ Ｂａｒｒｉｅｒ， ＨＦＢ）
子程序包处理。
２．４ 变密度地下水流及盐分运移模型
２．４．１ 数学模型 根据以上概念模型， 可建立描述研究区海水入侵过程的数学模型。 本文根据研究区
的水文地质条件及海水入侵的实际具体情形， 建立了龙口黄水河库区三维变密度地下水流和溶质运移
模拟模型， 采用 ＳＥＡＷＡＴ Ｖｅｒｓｉｏｎ ４（ＳＥＡＷＡＴ＿Ｖ４）计算软件对该模型进行了求解。 ＳＥＡＷＡＴ＿Ｖ４ 是由美国

地质勘探局（ＵＳＧＳ）开发的三维变密度地下水流和溶质运移模拟软件， 融合了地下水流模型 ＭＯＤＦＬＯＷ－
［２３］ 。 前者描述密度不断改变的混合咸淡水的流动； 后者描述地下水中
２０００ 和溶质运移模型 ＭＴ３ＤＭＳ
溶质（本研究中为盐分）的运移。 其中， 以等效淡水水头作为主因变量的描述三维变密度地下水流模型
的基本偏微分方程可写为：
æ é ｈ ρ－ρ ｚù ö æ é ｈ ρ－ρ ｚù ö æ é ｈ ρ－ρ ù ö ｈ ｆ ρ ｃ
ｆ
ｆ
ｆ
ｆ
ｆ
ｆ
ç ρＫ ê ＋ ú ÷ ＋ ç ρＫ ê ＋ ú ÷ ＋ ç ρＫ ê ＋ ú ÷ ＝ ρＳ ｆ ＋θ －ρ ｑ ｓ （１）
ｆｘ ê
ｓ
ú
ú
ｆｙ ê
ｆｚ ê
ú
ｘ è ë ｘ ρ ｘ û ø ｙ è ë ｙ ρ ｙ û ø ｚ è ë ｚ ρ û ø ｔ ｃ ｔ
ｆ
ｆ
ｆ
Ｈ（ｘ， ｙ， ｚ， ｔ） ＝ Ｈ （ｘ， ｙ， ｚ）， （ｘ， ｙ， ｚ）∈Ω （２）
ｔ ＝ ０ ０
Ｈ（ｘ， ｙ， ｚ， ｔ） ＝ Ｈ （ｘ， ｙ， ｚ， ｔ）， （ｘ， ｙ， ｚ）∈Γ ， ｔ≥０ （３）
Γ １ １ １
Ｈ
＝ Ｈ （ｘ， ｙ， ｚ， ｔ）， （ｘ， ｙ， ｚ）∈Γ ， ｔ≥０ （４）
ｎ ２ ２
Γ ２
式中： ｘ， ｙ， ｚ 分别为与主渗透方向一致的坐标轴； ρ 为淡水密度； ｈ 为等效淡水水头； Ｈ 为水头分
ｆ
ｆ
布函数； Ｈ 为初始时刻 ｔ ＝ ０ 时的地下水水头分布函数， 是已知函数； Ｋ 为等效淡水渗透系数； Ｓ 为等
０ ｆ ｆ
效淡水单位贮水系数； ρ 为流体密度； ρ 为源（汇）流体密度； ｑ 为单位体积源（汇）的流量； θ 为有效
ｓ
ｓ
孔隙度； ｃ 为溶质浓度； Ω 为研究区范围； 式（２）为初始条件； 式（３）为第一类边界条件， Γ 为第一类
１
边界； Ｈ （ｘ， ｙ， ｚ， ｔ）为边界 Γ 的已知函数； Γ 为第二类边界； Ｈ （ｘ， ｙ， ｚ， ｔ） 为边界 Γ 的已知
１ １ ２ ２ ２
函数。
在考虑溶质运移过程， 即盐分运移过程时， 选用化学性质上较为稳定的氯离子作为模拟组分。 在
海水入侵过程中， 氯离子的运移方式以对流和弥散为主， 暂不考虑运移过程中可能进行的化学反应。
描述研究区海水入侵区域氯离子运移行为的偏微分方程为：
ｃ ｃ ｃ （ｕ ｃ） （ｕ ｃ） （ｕ ｃ） ｃ ｑ ｓ
ｘ
ｙ
ｚ
∗
（Ｄ ｘｘ ）＋ （Ｄ ｙｙ ）＋ （Ｄ ｚｚ ）－ － － ＝ ＋ （ｃ－ｃ ） （５）
ｘ ｘ ｙ ｙ ｚ ｚ ｘ ｙ ｚ ｔ θ
Ｃ（ｘ， ｙ， ｚ， ｔ） ＝ Ｃ （ｘ， ｙ， ｚ）， （ｘ， ｙ， ｚ）∈Ω （６）
ｔ ＝ ０ ０
Ｃ（ｘ， ｙ， ｚ， ｔ） ＝ Ｃ （ｘ， ｙ， ｚ， ｔ）， （ｘ， ｙ， ｚ）∈Ｂ ， ｔ≥０ （７）
Ｂ １ １ １
ｎ ＝ ｑ （ｘ， ｙ， ｚ， ｔ）， （ｘ， ｙ， ｚ）∈Ｂ ， ｔ≥０ （８）
２ ２
－ＤｇｒａｄＣ· Ｂ ２
式中： Ｄ 为水动力弥散系数； μ 为实际流速； ｃ 为注水（抽水）的溶质浓度； Ｃ 为浓度分布函数； Ｃ 为
∗
０
初始时刻 ｔ ＝ ０ 时的浓度分布函数， 是已知函数； Ｂ 为一类边界； Ｃ （ｘ， ｙ， ｚ， ｔ）为一类边界 Ｂ 上的已
１ １ １
— １ ０ ９ —