一种地下水环境的建模及数值模拟方法 【技术领域】
本发明涉及一种数值模拟方法， 特别是关于一种地下水环境的建模及数值模拟方法。 背景技术 地下水环境突发污染事故需要快速掌握污染物在地下水环境中的迁移转化规律， 确定污染物扩散范围， 进而迅速实施有针对性的治理与防范。由此要求地下水环境数值模 型具有快速高效的建模与计算能力。
目前， 地下水环境的数值模拟方法主要有欧拉法、 拉格朗日法、 混合欧拉 - 拉格朗 日法等。基于这些方法构建了 MT3D、 HST3D、 WORM、 SWIFT、 USGS2D-MOC 等地下水环境数值模 拟软件。这些软件各有其优缺点， 但它们或因局限于饱和带或非饱和带单一含水层污染物 迁移模拟 ； 或因构建模型程序复杂， 难以适应快速建模的需要 ； 或因计算单元较多， 计算效 率低， 难以应用于地下水环境突发污染事故应急数值模拟等， 因此， 不能对地下水环境突发 污染事故进行全面、 高效、 准确的模拟计算。
发明内容
针对上述问题， 本发明的目的是提供一种地下水环境的建模及数值模拟方法， 该 方法可对地下水环境的突发污染事故进行全面、 高效、 准确的模拟计算。
为实现上述目的， 本发明采取以下技术方案 ： 一种地下水环境的建模及数值模拟 方法， 其包括以下步骤 ： 1) 前处理 ： 设置两个属性文件， 一是模拟区土壤属性文件， 包括研 究区边界、 X 与 Y 方向网格间距、 区域土壤渗透系数、 土壤有效孔隙度、 污染物在此种土壤中 的纵向和横向的弥散度、 土壤非饱和带扩散系数和溶质吸附系数 ； 二是地表高程与地下水 水位属性文件 ； 应用 IDL 软件的空间插值函数 KRIG2D() 或 TRI_SURF() 对地表高程、 地下水 位进行空间网格插值， 同时划分空间有限差分数值模拟网格， 得到两个相同的差分网格矩 阵： 地表高程矩阵 matrix1 和地下水位矩阵 matrix2 ； 应用 SIZE() 函数确定研究区网格属 性， 包括单元网格数、 单元网格行列数 ； 应用随机数值生成函数 FINDGEN() 生成随机数， 然 后将网格计算单元的节点坐标生成序列数组 ； 2) 二维水流模型计算二维水流模型计算包 括： 溶质非饱和带迁移模拟的地下水埋深矩阵计算以及饱和带二维地下水流场计算部分 ； 利用前处理得到的地表高程矩阵与地下水位矩阵相减得非饱和带厚度 h0， 即地下水位埋深 矩阵 matrixh0 ：
h0 ： matrixh0 ＝ matrix1-matrix2
该矩阵直接用于非饱和带溶质迁移模拟计算 ； 饱和带二维溶质迁移模型的水动力 条件是地下水流场的分布， 该流场由二维水流模型计算得到， 应用二维水流模型计算地下 水流场 ： ①计算单元节点流速 ： 利用达西定律 ：
其中， v 为流速 ； K 为渗透系数 ； I 为水力坡度 ； ne 为有效孔隙度 ； ②计算单元边界 面流速 ： 求得节点流速， 根据网格中相邻两个节点的流速， 取两个节点流速的平均值， 计算 得到网格单元边界面上流速分量的 IDL 矩阵 ； 3) 建立非饱和带一维溶质运移模型 ；
其中， t1 为溶质在非饱和带中的运移总时间 ； ΔH 为非饱和带厚度， ΔH ＝ h0 ； vz ＝ D/Δz 表示假设非饱和带厚度为 Δz ＝ 1m 时， 污染物垂向扩散速度， D 为垂向扩散系数， 量 纲为 L2/T ； 4) 建立饱和带溶质迁移二维数值模型 ：
其中， R 为等温吸附系数， C 为溶解浓度， 单位为 ML-3， t1 为饱和带中溶质迁移数值 2 -1 模拟的时间， Dij 为弥散系数张量， 单位为 L T ， x、 y 为计算距离 L， vx 和 vy 为地下水流速， -1 3 -1 单位为 LT ， ， qs 为源 / 汇处单位体积含水层流量 (， 单位为 M T ， θ 为含水层孔隙度， Cs 为 -3 -1 源 / 汇项浓度， 单位为 ML ， λ1 为溶解相的反应速率常数， 单位为 T ， λ2 为吸附相的反应 速率常数， 单位为 T-1， ρb 为空隙介质体积密度， 单位为 ML-3， 为吸附浓度， 单位为 MM-1 ； 5) 将公式 (2) 中的各项式进行近似表示 ： 在规则间距的节点中心网格中， 公式 (2) 右边第一项 在单元 (i， j) 近似为 ：
公式 (3) 表示 x 方向上， 由浓度梯度引起从 x 方向上进入单元 (i， j) 的净弥散通 量； 公式 (2) 右边第二项近似为 ：
公式 (4) 表示由 y 方向的浓度梯度引起的从 x 方向进入单元 (i， j) 的净弥散通 量； 同理， 得到由 y 方向的浓度梯度引起的从 y 方向上进入单元 (i， j) 的净弥散通量， 即公 式 (2) 右边的第三项为 ：
由 x 方向的浓度梯度引起的从 y 方向进入单元 (i， j) 的净弥散通量， 即公式 (2) 右边的第四项， 写为 ：
公式 (2) 右边第五项为从 x 方向进入单元 (i， j) 的净对流通量， 近似为 ：第六项为从 y 方向进入单元 (i， j) 的净对流通量， 近似为 ：其中， α 为空间加权因子 ； 中心加权， α 等于 0.5 ； 采用上风加权方案 ：公式 (2) 右边第七项是由源 / 汇引起的流入或流出单元 (i， j) 净质量通量 ； 公式 (2) 右边第八项和第九项是单元 (i， j) 内一级化学反应引起的质量损失或增加量 ； 吸附浓 度 表示为溶解浓度的函数， 具体取决于所用的吸附等温线 ；
公式 (2) 的左边是储存于单元 (i， j) 的质量变化率， 近似为 ：
其中延迟因子 R 由吸附等温线确定 ； 6) 对以上各分量差分计算过程详细实现过程 如下 ： ①弥散系数 D ： 由边界流速分量及公式 (10)， 利用 IDL 矩阵运算计算弥散系数分量 ；
其中 ： D* 为有效分子扩散系数 ；②弥散项 ： 设置初始溶质浓度场， 即确定污染事件发生位置及该物质的饱和溶解度作为模型的初始浓度分布矩阵 C ； 由弥散项计算公式 (3)、 公式 (4)、 公式 (5) 及公式 (6)， 结合弥散系数计算结果矩阵， 利用 IDL 矩阵运算计算各弥散项 ； ③对流项 ： 由对流项计算 公式 (7)、 公式 (8)， 结合单元格边界流速分量、 土壤弥散度属性， 利用 IDL 矩阵运算计算各 对流项 ； ④公式 (2) 的第七项， 源汇项 ： 由于快速处理模型、 源汇项设置为 0 ； ⑤公式 (2) 的 第八、 第九项， 化学反应项 ： 根据各化学物质的等温吸附曲线获取 ； ⑥单元节点浓度 ： 由公 式 (9) 及各个分量矩阵计算结果， 利用 IDL 矩阵运算计算单元节点后一时刻浓度场 ； 7) 返 回步骤 5) 中的计算过程， 即实现一个时间段后， 区域流场及溶质浓度场的计算 ； 在步骤 5) 前加一个总时段的循环控制语句， 即实现各时段区域流场及溶质浓度场的计算 ； 8) 结果输 出： 经过上述步骤实时输出地表高程矩阵 matrix1、 地下水位矩阵 matrix2、 地下水埋深矩 阵 matrixh0、 地下水流场、 溶质浓度场和时间项。
所述步骤 4) 中， 根据中心加权方案， 从节点浓度获得界面浓度 ：
本发明由于采取以上技术方案， 其具有以下优点 ： 1、 本发明建立的一维模型和二 维模型的结构简单， 模型核心计算全部采用矩阵计算， 矩阵公式与算法公式完全一致， 易于 理解与建模。2、 本发明由于充分运用了 IDL 平台提供的函数库， 对于空间插值、 矩阵计算等 原本在 Fortran 等平台下十分繁琐的工作， 本模型多数仅需一行代码即可完成， 全部模型 仅有约 180 行代码， 模型代码精炼。3、 本发明由于核心处理采用了基于 IDL 矩阵及函数计 算， 模型计算效率大为提高。4、 本发明针对地下水环境突发污染事件， 选定模拟区域后， 可 快速插值生成指定大小的计算网格 ； 并根据污染区域周边少数观测井的地表高程、 潜水水 位、 地表水水位， 快速插值生成近似于实际条件的地形、 水位及水位埋深等 ； 通过快速建立 的模型可指定污染物位置、 浓度等。 5、 本发明采用有限差分方法计算单元溶质的迁移扩散， 完全摆脱了以往的 Fortran 语言或 C 语言的数组循环计算， 取而代之的是 IDL 平台下的矩阵运算方法， 计算简单高效。本发明基于 IDL 平台， 可高效、 快速、 简便地模拟地下水环境， 因此， 可广泛用于地下水环境突发污染事件的数值模拟、 预测、 应急预测等过程中。 附图说明
图 1 是本发明流程示意图
图 2 是本发明计算由 x 以及 y 方向的浓度梯度引起的 x 方向上的弥散通量示意图
图 3 是本发明快速处理模块地表高程插值结果示意图
图 4 是本发明快速处理模块地下水位插值结果示意图
图 5 是本发明快速处理模块潜水埋深计算结果示意图
图 6 是本发明实施例中地下水环境应急快速预警模型模拟结果示意图
图 7 是本发明实施例中当 t ＝ 1000 天时区域流场与污染场叠加分布示意图 具体实施方式
下面结合附图和实施例对本发明进行详细的描述。
本发明是一种耦合土壤非饱和带、 饱和带溶质迁移过程的地下水环境数值模型的 快速构建与准确模拟方法， 其基于以下思想 ： 由于污染事件通常发生于地表， 污染物首先经过含水层上部非饱和带垂向运动， 然后进入含水层饱和带 ； 本发明对在非饱和带含水层中 的溶质迁移采用垂向一维模型， 而在饱和含水层中， 由于溶质垂向迁移扩散速率往往较横 向迁移扩散速率小一个数量级， 因此， 基于快速模拟计算忽略饱和含水层的垂向运动， 采用 了水平二维数值模拟方法。
本发明包括以下步骤， 如图 1 所示 ：
1) 前处理 ：
设置两个属性文件， 一是模拟区土壤属性文件 (Soilpro 文件 )， 包括研究区边界 ( 模拟区范围 )、 X 与 Y 方向网格间距、 区域土壤渗透系数、 土壤有效孔隙度、 污染物在此种 土壤中的纵向和横向的弥散度、 土壤非饱和带扩散系数和溶质吸附系数， 如表 1 所示 ； 二是 地表高程与地下水水位属性文件 (Wellpro 文件 )， 如表 2 所示。
应用 IDL 软件的空间插值函数 KRIG2D() 或 TRI_SURF() 对地表高程、 地下水位进 行空间网格插值， 同时划分空间有限差分数值模拟网格， 得到两个相同的差分网格矩阵 ： 地 表高程矩阵 (matrix1) 和地下水位矩阵 (matrix2)。应用 SIZE() 函数确定研究区网格属 性， 包括单元网格数、 单元网格行列数等。应用随机数值生成函数 FINDGEN() 生成随机数， 然后将网格计算单元的节点坐标生成序列数组。
表 1Soilpro 的文件格式
表 2Wellpro 的文件格式0 经度 X 1 维度 Y 2 地表高程 ( 米 ) 3 地下水位 ( 米 )2) 二维水流模型计算
二维水流模型计算包括 ： 溶质非饱和带迁移模拟的地下水埋深矩阵计算以及饱和 带二维地下水流场计算部分。
利用前处理得到的地表高程矩阵与地下水位矩阵相减可得非饱和带厚度 h0， 即地 下水位埋深矩阵 matrixh0 ：
h0 ： matrixh0 ＝ matrix1-matrix2
该矩阵可直接用于非饱和带溶质迁移模拟计算。
饱和带二维溶质迁移模型的水动力条件是地下水流场的分布， 该流场可由二维水 流模型计算得到， 应用二维水流模型计算地下水流场 ：
①计算单元节点流速
利用达西定律 ：
其中， v 为流速 ； K 为渗透系数 ； I 为水力坡度 ； ne 为有效孔隙度。 根据达西定律和表 1， 利用 IDL 矩阵计算节点流速代码如下 ：②计算单元边界面流速
求得节点流速， 根据网格中相邻两个节点的流速， 取两个节点流速的平均值， 计算 得到网格单元边界面上流速分量的 IDL 矩阵， 计算代码如下 ：
3) 非饱和带一维溶质运移模型计算
非饱和带计算采用简化一维溶质运移模型， 仅计算溶质由地表运移至含水层的时 间， 计算公式如下 ：
其中， t1 为溶质在非饱和带中的运移总时间 ； ΔH ＝ h0 为非饱和带厚度 ； vz ＝ D/ Δz 表示假设非饱和带厚度为 Δz ＝ 1m 时， 污染物垂向扩散速度， D 为垂向扩散系数， 量纲 2 为 (L /T)。
4) 饱和带溶质迁移二维数值模型计算
采用有限差分法的欧拉方法， 计算地下水溶质迁移扩散过程。模型饱和含水层溶 质对流 - 弥散迁移采用二维对流 - 弥散方程如下 ：
其中， R 为等温吸附系数， C 为溶解浓度 (ML-3)， t1 为饱和带中溶质迁移数值模拟 2 -1 的时间， Dij 为弥散系数张量 (L T )， x、 y 为计算距离 (L)， vx 和 vy 为地下水流速 (L T-1)， ，qs 为源 / 汇处单位体积含水层流量 (M3T-1)， θ 为含水层孔隙度， Cs 为源 / 汇项浓度 (ML-3)， λ1 为溶解相的反应速率常数 (T-1)， λ2 为吸附相的反应速率常数 (T-1)， ρb 为空隙介质体 积密度 (ML-3)， 为吸附浓度 (MM-1)。该模型充分考虑了弥散、 迁移、 源汇项、 化学吸附等影 响。
5) 将公式 (2) 中的各项式进行近似表示， 以便采用有限差分法编程实现。 在规则间距的节点中心网格中， 公式 (2) 右边第一项在单元 (i， j) 近似为 ：如图 2(a) 所示， 公式 (3) 表示 x 方向上， 由浓度梯度引起从 x 方向上进入单元 (i， j) 的净弥散通量。
公式 (2) 右边第二项可近似为 ：
如图 2(b) 所示， 公式 (4) 表示由 y 方向的浓度梯度引起的从 x 方向进入单元 (i， j) 的净弥散通量。其中用到中心加权方案， 从节点浓度获得界面浓度， 图中以 标记出 ：
同理， 可以得到由 y 方向的浓度梯度引起的从 y 方向上进入单元 (i， j) 的净弥散 通量， 即公式 (2) 右边的第三项为 ：
由 x 方向的浓度梯度引起的从 y 方向进入单元 (i， j) 的净弥散通量， 即公式 (2) 右边的第四项， 可写为 ：
公式 (3)、 公式 (4)、 公式 (6) 及公式 (7) 中的弥散系数分量的计算， 若给定纵向弥 散度 αL 与横向弥散度 αT， 各弥散系数计算公式如下 ：
其中 ： D* 为有效分子扩散系数。公式 (2) 右边第五项为从 x 方向进入单元 (i， j) 的净对流通量， 近似为 ：第六项为从 y 方向进入单元 (i， j) 的净对流通量， 可近似为 ：其中 α 为空间加权因子。中心加权， α 等于 0.5 ； 采用上风加权方案 ：公式 (2) 右边第七项是由源 / 汇引起的流入或流出单元 (i， j) 净质量通量， 第八 项和第九项是单元 (i， j) 内一级化学反应引起的质量损失或增加量。吸附浓度 可以表示 为溶解浓度的函数， 具体取决于所用的吸附等温线。
公式 (2) 的左边是储存于单元 (i， j) 的质量变化率， 可以近似为 ：
其中延迟因子 R 由吸附等温线确定。 6) 对以上各分量差分计算过程详细实现过程如下 ： ①弥散系数 D ： 由边界流速分量及公式 (8)， 利用 IDL 矩阵运算可以计算弥散系数分量 ； ②弥散项 ： 设置初始溶质浓度场， 即确定污染事件发生位置及该物质的饱和溶解 度作为模型的初始浓度分布矩阵 C。 由弥散项计算公式 (3)、 公式 (4)、 公式 (6) 及公式 (7)， 结合弥散系数计算结果矩阵， 利用 IDL 矩阵运算可以计算各弥散项 ；
③对流项 ： 由对流项计算公式 (10)、 公式 (11)， 结合单元格边界流速分量、 土壤弥 散度属性， 利用 IDL 矩阵运算可以计算各对流项 ；
④公式 (2) 的第七项， 源汇项 ： 由于快速处理模型， 源汇项可以暂时设置为 0 ；
⑤公式 (2) 的第八、 第九项， 化学反应项 ： 根据各化学物质的等温吸附曲线获取 ；
⑥单元节点浓度 ： 由公式 (12) 及各个分量矩阵计算结果， 利用 IDL 矩阵运算可以 计算单元节点后一时刻浓度场 ；
计算单元节点后一时刻浓度场的 IDL 矩阵计算代码如下 ：
cinit[1:nDG[0]-1， 1:nDG[1]-1] ＝ cinit[1:nDG[0]-1， 1:nDG[1]-1]+$
(deltT/RR)*(Dcxx[0:*， 0:*]+Dcxy[0:*， 0:*]+Dcyy[0:*， 0:*]+$
Dcyx[0:*， 0:*]+Vqcx[0:*， 0:*]+Vqcy[0:*， 0:*]+Lambda)
其中， cinit[] 为浓度场矩阵， deltT 为时间步长， RR 为延迟因子， Dcij[] 为弥散 项， Vqcx[] 为对流项， Lambda 为等温吸附项。
7) 然后返回步骤 5) 中的计算过程， 即可实现一个时间段后， 区域流场及溶质浓度 场的计算 ； 在步骤 5) 前加一个总时段的循环控制语句， 即可实现各时段区域流场及溶质浓 度场的计算。
8) 结果输出 ： 经过上述步骤可以实时输出地表高程矩阵 matrix1、 地下水位矩阵 matrix2、 地下水埋深矩阵 matrixh0、 地下水流场、 溶质浓度场和时间项 ( 时间段 )。
上述步骤 8) 输出的结果以及二维水流模型结果可与地表水环境模型接口连接， 用于地表水环境， 数值模拟河道侧向入流地下水的溶质边界条件， 且输出结果可根据地表 水环境模型输入文件的格式需要设计输出相应的文件格式。
下面列举一具体实施例 ：
1、 水流模拟
假设任意一个地方发生污染事件， 利用前处理模块可以快速得到模型研究区域边 界， 设定模型剖分网格精度、 土壤属性等， 如表 3 所示， 为土壤属性即 Soilpro 的文件。根据 研究区提供的 5 眼地下水位观测井的资料， 如表 4 所示， 为地下水位和地表高程即 Wellpro 的文件， 可以插值得到该区域地表高程， 如图 3 所示。地下水水位插值结果， 如图 4 所示 ； 地 下水埋深插值结果， 如图 5 所示。 图 3 ～ 5 中的 (a) 图为采用 Kringing 方法的差值结果图， (b) 图所示为采用五次样条光滑差值方法的差值结果图。
表 3Soilpro 的文件
表 4Wellpro 的文件13101908100 A CN 101908105
说0 经度 X 20229711 20230350 20237854 20236626 20243586 1 维度 Y 3598113 3603662 3603889 3597749 3595293 2明书3 地下水位 ( 米 ) 28 30.5 30.1 27.5 27.210/10 页地表高程 ( 米 ) 30.5 33.4 32.3 30.1 29.72、 溶质迁移模拟
假设污染发生位置坐标为 (20238486， 3599112)， 假设该污染物在含水层中饱和溶 解度为 10000mg/L。通过模拟计算， 可得到以下结果 ：
(1) 污染物由地表运移通过非饱和带时间为 32.11 天。 该处地下水埋深为 6.422m。
(2) 污染物进入含水层 5 天后向东北方向最远扩散距离为 140m( 扩散标准为区域 污染物浓度达到 0.01mg/L) ； 污染物进入含水层 10 天后向东北方向最远扩散距离为 141m ； 污染物进入含水层 30 天后向东北方向最远扩散距离为 214m ； 污染物进入含水层 100 天后 向东北方向最远扩散距离为 278m ； 污染物进入含水层 500 天后向东北方向最远扩散距离为 390m ； 污染物进入含水层 1000 天后向东北方向最远扩散距离为 450m。各时间节点污染物 浓度分布图， 如图 6 所示。
(3) 通过地下水流场与第 1000 天污染物浓度场叠加图， 如图 7 所示， 可以看出， 模 型模拟污染物迁移扩散趋势完全符合流场变化规律。
上述实例的关键技术在于 ：
实例研究中， 研究区共计剖分 51506 个计算单元， 计算时间步长 1 天， 时间总长度 为 1000 天， 模型输出为 1000 天的各节点地表高程、 水位、 水位埋深、 溶质浓度。计算机硬件 配置为 Intel(R)Core(TM)2 Quad Q6600 2.40GHz CPU， 内存 4GB(PC2-5300 DDR2 666MHz)， 硬盘为希捷 ST3500820AS(498GB， 7200 转 / 分 )。操作系统为 Windows XP(32bit/SP2)。模 型计算耗时 50.18 秒。
上述各实施例仅用于说明本发明， 其中各部件的结构、 连接方式等都是可以有所 变化的， 凡是在本发明技术方案的基础上进行的等同变换和改进， 均不应排除在本发明的 保护范围之外。