在本节中，运用Processing Modflow软件，对研究区内第四系地下水系统地下水流场变化情况进行数值模拟计算，研究第四系地下水系统流场的演变情况，本次模拟总面积与研究区面积一致，为25km2。通过数值模拟，更加细致地刻画了研究区第四系地下水系统中六价铬污染晕的变化情况，进而为研究抽水处理技术的可行性奠定基础。
11.3.2.1 水文地质概念模型
(1)含水层和弱透水层概化
水文地质概念模型是研究区水文地质条件的综合和概化，是建立数学模型的基础，根据研究区水文地质条件及钻孔资料分析，研究区内第四系地下水系统主要由多层交互的砂、粉砂、粉土与黏土层构成，含水层组岩相变化与分布规律明显。从地下水流动系统的观点看，本区第四系地下水系统的补给、径流与排泄条件清楚;地下水除了在水平方向上发生交换外，垂向上的交换是必然的。
根据水文地质特征，结合含水层组的划分，通过一定的钻孔资料分析，将第四系地下水系统概化为自上向下的含水层(第Ⅱ含水组)→弱透水(隔水层)→含水层(第Ⅲ含水组)结构系统(图11.24)。含水层岩性以砂性土为主，含有若干相对较薄的不连续的黏性土夹层或透镜体，但每个含水层仍然具有较强的统一的水力联系。弱透水层主要是由厚度较大的连续的黏性土组成，构成承压含水层的顶板和底板。相邻含水层在水头差作用下，通过其发生越流，含水层与弱透水层构成多层结构的含水系统。
图11.24 研究区第四系地下水系统水文地质概念模型
(2)初始条件及边界条件的概化
据研究区内第四系地下水系统地下水流特征及地层结构分析，将侧向边界确定为流量边界;其量根据第四系地下水系统地下水流场、地层结构和周边地下水动态观测井的情况来确定。第四系地下水系统水文地质概念模型的第一层顶部为模型的上边界，通过该边界，与外界发生垂向水量交换，如大气降水入渗补给、河流与灌溉入渗补给、潜水的蒸发及植物的蒸腾作用等。底部以第Ⅲ含水组底部为边界，由于第Ⅲ含水层组底部为第Ⅳ含水层组，第Ⅳ含水层组赋存于新生界系;该组岩性在上部有较厚的泥岩和泥砾岩层，成为隔断第Ⅲ含水组与第Ⅳ含水层组水力联系的隔水层。第四系地下水系统水文地质概念模型第三层和第一层之间通过第二层越流进行水量交换。
11.3.2.2 数学模型的建立与求解
(1)数学模型的建立
分析研究区第四系地下水系统水文地质结构模型，在确立了地下水侧向、垂向边界后，根据已掌握地下水动态特征，确定研究区为非均质各向同性非稳定流地下水系统，其数学模型[62]如下:
潜水含水层数学模型为
变环境条件下的水资源保护与可持续利用研究
变环境条件下的水资源保护与可持续利用研究
变环境条件下的水资源保护与可持续利用研究
变环境条件下的水资源保护与可持续利用研究
承压含水层数学模型为
变环境条件下的水资源保护与可持续利用研究
变环境条件下的水资源保护与可持续利用研究
变环境条件下的水资源保护与可持续利用研究
变环境条件下的水资源保护与可持续利用研究
式中:H1、H2为潜水和承压的水头;K1为潜水含水层渗透系数;T1为承压含水层导水系数;K'、M'为潜水含水层于承压含水层之间弱透水层或岩性“天窗”处介质的垂向渗透系数及厚度;B为潜水含水层底板标高;W为源汇项;ε为承压水含水层的开采强度;μe、μd为潜水含水层的重力给水度和承压水含水层的弹性给水度;t为时间;qn为单位流量在垂直Г上的分量;h0为初始水位;Ω1、Ω2为潜水含水层和承压水含水层计算区范围;Г1为已知水位边界(一类边界);Г2为已知流量边界(二类边界)。
(2)数学模型的求解
本次研究采用地下水数值模拟软件Processing Modflow，利用美国地质调查局所开发的三维有限差分地下水流模型Modflow和IBS软件包进行模拟求解。Processing Modflow是由 Wen-Hsing Chiang和 Wolfgang Kinzelbach等在以 Modflow为基础研制开发的模拟地下水运动和溶质运移的计算软件，能够模拟由于抽水处理技术抽取地下水引起的含水层的压缩量，由Interbed-Storage Package(IBS)来实现。
其求解方法是在计算区域内采用矩形剖分和线性插值，应用有限差分法将上述数学模型离散为有限单元方程组，然后求解。
(3)空间离散
整个模拟研究区面积为25km2，据水文地质概念模型范围、第四系地下水系统含水层结构特征及第四系地下水系统地下水流动特征，将研究区剖分为69行104列，共剖分7176个单元格(图11.25)。
图11.25 研究区第四系地下水系统网格剖分立体示意图
(4)时间离散
地下水渗流过程是个缓慢过程，往往需要较长时间序列的观测资料，进行综合分析。研究区第四系地下水系统地下水流场计算利用资料最为齐全的2005～2009年水位观测资料进行模拟，模拟期为2005年1月到2009年12月末，将整个模拟期划分为59个应力期，计算时间步长为1个月(30d)，预测期为2010年12月～2020年12月。
11.3.2.3 汇源项的确定
(1)降水入渗补给量
研究区第四系地下水系统最主要的补给来自大气降水。降雨量主要通过气象站的观测数据获得，与相应的降雨入渗系数(图11.26)进行属性叠加计算，得到降雨入渗量。除此之外，该地区的补给还包括河流入渗量，农业灌溉入渗量和由第四系地下水系统流场及达西定律计算得到侧向流量等(表11.8)。计算公式如下:
Q降=αkPF (11.13)
式中:Q降为降水入渗补给总量，m3/a;α为降水入渗系数，根据计算区表层土性质和多年水位埋深综合选取;k为有效入渗面积调整系数;P为计算区多年平均降雨量，m/a;F为计算区面积，m2。
(2)越流补给量
研究区内第四系地下水系统第Ⅱ含水组和第Ⅲ含水组之间存在一定的水位差，两个含水组之间隔水层相对较薄，累积厚度为6～8m;再加上由于研究区内大量施工供水井，人为沟通了两个含水组，存在第四系地下水系统第Ⅱ含水组对第四系地下水系统第Ⅲ含水组的越流补给，越流补给量计算公式如下:
变环境条件下的水资源保护与可持续利用研究
式中:Qy为相邻含水层的越流补给量，m3/d;Kμ、Ke为开采层上、下部弱透水层垂直渗透系数，m3/d;Mμ、Me为开采层上、下部弱透水层厚度，m;Fμ、Fe为开采层上、下部越流面积，m;Hμ、He为开采层的越流层水位，m;h为开采层的水位或开采漏斗的平均水位，m。
图11.26 研究区第四系地下水系统降水入渗系数分区图
表11.8 研究区第四系地下水系统补给量和排泄量计算表
(3)排泄量
研究区内第四系地下水系统地下水排泄方式主要是人工开采。人工开采主要包括B钢厂水源地，A钢厂水源地和城区第一水源地的开采;另外，研究区内其他自备井分布比较均匀，故将其地下水开采量总量进行平均，然后均匀分布于全区进行计算。
通常认为水位埋深大于4m的地区潜水蒸发量较小，而研究区内第四系地下水系统潜水分布于地表5m以下，可认为研究区蒸发量较小(表11.8)。
Q蒸=ckεF (11.15)
式中:Q蒸为蒸发排泄总量，m3/a;c为蒸发系数，根据计算区表层土性质和多年水位埋深综合选取;k为有效蒸发面积调整系数;ε为计算区多年平均蒸发量，m/a;F为计算区面积，m2。
11.3.2.4 模型的识别与验证
(1)初始流场的确定
研究区内第四系地下水系统第Ⅱ含水组长观孔共5个，第Ⅲ含水组长观孔共6个(图11.27)。研究区第四系地下水系统地下水流场模拟选择2005年1月研究区第四系地下水系统地下水等水位线图作为初始流场(图11.28和图11.29)。
图11.27 研究区内第四系地下水系统长观孔分布图
图11.28 研究区2005年1月第四系地下水系统第Ⅱ含水组地下水等水位线图
图11.29 研究区2005年1月第四系地下水系统第Ⅲ含水组地下水等水位线图
(2)水文地质参数的确定
研究区第四系地下水系统地下水流模拟模型水文地质参数包括第四系地下水系统第Ⅱ和第Ⅲ含水层组及弱透水层的渗透系数、重力给水度及释水系数。据收集的所在区“地下水资源开发区划报告”导水系数分区图和多组抽水试验所取得的各井孔第四系地下水系统第Ⅱ和第Ⅲ含水层组的水文地质参数作为数值模型计算的初步参考，将研究区第四系地下水系统第Ⅱ和第Ⅲ含水层组的水文地质参数分为以下几个量值区(表11.9和表11.10，图11.30和图11.31)。
表11.9 研究区第四系地下水系统第Ⅱ含水层组渗透系数分区取值表
表11.10 研究区第四系地下水系统第Ⅲ含水层组渗透系数分区取值表
图11.30 研究区第四系地下水系统第Ⅱ含水层组水文地质参数分区示意图
图11.31 研究区第四系地下水系统第Ⅲ含水层组水文地质参数分区示意图
(3)模型验证结果
从2009年11月第四系地下水系统第Ⅱ含水组(5个长观孔)和第Ⅲ含水组(6个长观孔)等水位线拟合曲线图(图11.32和图11.33)可以看出，研究区第四系地下水系统地下水等水位线与实际较为吻合，计算的第四系地下水系统地下水流场变化与实际地下水流场的变化在数量上基本一致。
图11.32 研究区2009年11月第四系地下水系统第Ⅱ含水组等水位线拟合曲线图
图11.33 研究区2009年11月第四系地下水系统第Ⅲ含水组等水位线拟合曲线图
图11.34 研究区第四系第Ⅱ含水组观05地下水水位观测值与计算值拟合曲线图
选取研究区第四系地下水系统第Ⅱ含水组具有典型代表性的观05、观08和第Ⅲ含水组具有典型代表性的TC004、TC005观测孔进行拟合，研究区第四系地下水系统第Ⅱ含水组水位动态观测值与计算值拟合曲线图(图11.34和图11.35)曲线变化相对较为平缓，实测值曲线各年波动相对较大;第Ⅲ含水组水位动态观测值与计算值拟合曲线图(图11.36和图11.37)曲线变化相对较为平稳，实测值曲线各年波动相对较大。动态结果表明第四系地下水系统第Ⅱ含水组和第Ⅲ含水组地下水位动态观测值与计算值变化趋势一致;经过不断调参，所建立的第四系地下水系统地下水流数学模型对水文地质条件的概化、边界条件的选择以及参数的确定等都与研究区的实际情况符合较好。最终确定研究区第四系地下水系统水文地质参数(表11.11、表11.12)，可以作为该模型地下水水位及流场预报的参数。
图11.35 研究区第四系第Ⅱ含水组观08地下水水位观测值与计算值拟合曲线图
图11.36 研究区第四系第Ⅲ含水组TC004地下水水位观测值与计算值拟合曲线图
图11.37 研究区第四系第Ⅲ含水组TC005地下水水位观测值与计算值拟合曲线图
表11.11 研究区第四系地下水系统第Ⅱ含水层组渗透系数分区取值表
表11.12 研究区第四系地下水系统第Ⅲ含水层组渗透系数分区取值表
11.3.2.5 第四系地下水系统地下水流场预测
预测时间为2010年11月至2020年12月;其中2011年6月以后，A钢厂和B钢厂水源地全部关闭。根据第四系地下水系统地下水流模型进行预测，并分别给出两大水源地关停后2011年年底、2012年年底、2015年年底和2020年年底第四系地下水系统第Ⅱ含水组和第Ⅲ含水组地下水流场预测图(图11.38～图11.45)。
图11.38 研究区2011年12月第四系地下水系统第Ⅱ含水组地下水流场预测图
图11.39 研究区2012年12月第四系地下水系统第Ⅱ含水组地下水流场预测图
图11.40 研究区2015年12月第四系地下水系统第Ⅱ含水组地下水流场预测图
图11.41 研究区2020年12月第四系地下水系统第Ⅱ含水组地下水流场预测图
图11.42 研究区2011年12月第四系地下水系统第Ⅲ含水组地下水流场预测图
图11.43 研究区2012年12月第四系地下水系统第Ⅲ含水组地下水流场预测图
图11.44 研究区2015年12月第四系地下水系统第Ⅲ含水组地下水流场预测图
图11.45 研究区2020年12月第四系地下水系统第Ⅲ含水组地下水流场预测图
另外，分别将现状年(2010年)第四系地下水系统第Ⅱ含水组和第Ⅲ含水组的六价铬污染晕与预测的地下水流场形态进行叠加，以便分析六价铬污染晕的发展趋势。
(1)第四系地下水系统Ⅱ含水组地下水流场预测
对比研究分析以上第四系地下水系统第Ⅱ含水组地下水流场预测图可得:
1)水位整体上升。关闭两个钢厂水源地后，研究区第四系地下水系统Ⅱ含水组地下水位逐年上升。2011年年底，地下水水位为-19.1～-7.8m;2012年年底，地下水水位为-13.8～-6.9m;2015年年底，地下水水位为-7.9～-3.7m，到了2020年年底，地下水水位达到1.4～3.5m。研究区内第四系地下水系统两大水源地抽水井关闭，极大地减少了第Ⅱ含水组的人工排泄量，引起第Ⅱ含水组地下水位逐年上升。
2)地下水位降落漏斗中心转移，地下水流动速度变缓。第四系地下水系统Ⅱ含水组地下水位降落漏斗中心向西转移至城区第一水源地。2011年年底的地下水位降落漏斗中心水力梯度为46‰，2012年年底的地下水位降落漏斗中心水力梯度为1.5‰，2015年年底的地下水位降落漏斗中心水力梯度为0.5‰，到2020年年的地下水位降落漏斗中心水力梯度底仅为0.08‰。以上说明漏斗中心水力坡度(地下水流动交换速度)逐年减小，至2020年，地下水位降落漏斗将基本消失。
3)污染晕呈扩展趋势。两钢厂水源地关闭以后，研究区第四系地下水系统Ⅱ含水组地下水流场的变化主要表现在地下水水位上升和水流速度减缓。另外，研究区区域地下水水流方向也会发生一定的变化，但污染晕范围内地下水流向变化不大。结合现状污染晕形态和地下水流方向的变化可知，2015年以后，第四系地下水系统Ⅱ含水组污染晕会继续向西扩展，即向丰南第一水源地方向扩展。
(2)第四系地下水系统Ⅲ含水组地下水流场预测
对比研究分析由以上第四系地下水系统第Ⅲ含水组地下水流场预测图可得:
1)水位整体上升。关闭两钢厂水源地抽水井，导致研究区内第四系地下系统第Ⅲ含水组排泄量大幅下降，以致地下水水位逐年上升;2012年年底，地下水水位为-55.6～-27.3m;2013年年底，地下水水位为-50.1～-24.5m;2015年年底，地下水水位为-45.1～-21.7m;到2020年年底，地下水水位达到-43.2～20.3m，表明研究区第Ⅲ含水组地下水水位整体呈现较大幅度的上升。
2)地下水水位漏斗中心转移，地下水流动速度变缓。第四系地下系统第Ⅲ含水组地下水水位降落漏斗中心从A钢厂北厂院内逐渐转移至城区第一水源地，说明在关停水源地后，第Ⅲ含水组地下水向城区第一水源地流向汇集。2011年年底，地下水水位降落漏斗中心水力梯度为62.5‰;2012年底，地下水水位降落漏斗水力梯度为32.7‰;2015年底，地下水水位降落漏斗水力梯度为12.7‰;到2020年年底，地下水水位降落漏斗水力梯度仅为6.0‰。说明漏斗中心水力坡度(地下水流动交换速度)逐年减小，漏斗中心范围逐渐缩小;但至2020年，地下水位降落漏斗不会消失。
3)污染晕呈扩展趋势。钢厂水源地关闭后，研究区第四系地下水系统第Ⅲ含水组地下水流场的变化主要表现在水位上升和区域水流方向发生一定的变化，水流速度也会减缓。结合现状年(2010年)污染晕形态和地下水流方向的变化可知，2011年以后，第Ⅲ含水组污染晕会继续向西扩展，即向城区第一水源地方向扩展，其扩展速度比第Ⅱ含水组污染晕要快。