数值法水量计算的步骤

数值法水量计算的正式进行，是在水文地质勘察外业工作完成之后，计算所需的各种原始数据按照数值法水量计算的要求已全部整理、分析、计算完毕，只等调用；计算所需的任何一种数据缺失，都会导致水量计算无法进行。

一般而言，数值法水量计算需要经过如下几个重要的步骤。

1.概化水文地质模型

包括对计算区域的圈定、计算域边界形状及性质的确定、含水层的划分、上下层水力联系的确定、参数分区等。

其中，计算区域的圈定，主要考虑计算任务的需要，兼顾边界条件处理的方便。对于数百、数千、数万甚至数十万平方千米的大区域，多在勘察区边界附近分段选定出方便概化处理的计算区边界即可，计算区边界与勘察区边界两者一般不会相差非常悬殊。然而，对于面积仅有数平方千米、或者最多数十平方千米的供水水源地，计算区边界圈定与水源地范围常常悬殊很大，并且必须比水源地范围要大。根据情况不同，常会作这样的处理：

（1）当近距离内没有已建成的相邻水源地时，计算区域的圈定主要考虑本水源地开采方案布设的方便、边界条件概化的方便。

（2）当近距离内有已建成的相邻水源地、新建水源地开采有可能影响到相邻水源地时，计算区域的圈定除了考虑本水源地开采方案布设的方便、边界条件概化的方便之外，还应当把已建成的相邻水源地也纳入计算域内，计算本水源地布设的开采方案时，把相邻水源地现状的开采方案也纳入计算中，以论证计算本水源地布设的不同开采方案对相邻水源地的影响。

（3）当计算区的某一侧是无限延伸时，也就是说开采方案引起的流场变化远远达不到这一侧的自然边界时，如果也不具备一条已知水头变化规律的“已知变水头的边界”时，可以把这一侧的模型边界向外推至开采方案引不起流场变化的数公里之外，作为已知的定水头边界处理。

2.选用合适的数学模型来描述

是潜水还是承压水还是承压转无压水，或者是二者、三者都有；单层还是多层含水层；各向同性还是各向异性；二维流还是准三维流还是三维流；稳定流还是非稳定流；以及初始条件、边界条件的数学描述。

3.单元剖分

在参数分区的基础上进行单元剖分，单元剖分时应当注意：①剖分的单元不能一个单元骑跨在两个参数区上；②用作拟合孔、预测孔的观测孔要剖分在节点上。

4.编程

自己编制计算程序，或者购买适合的计算软件。如果是自己编制计算程序，则要经过调试程序、检验程序，确认程序无误后，再投入正式使用。如果是购买适合的计算软件，则要经过学习、使用软件的培训。

5.计算时段剖分

先确定模型识别、模型检验所涉及起始日期、时间，结束日期、时间，然后，按照一定的递推关系进行计算时段的剖分。

6.初始流场、末刻流场的模拟

在统测水位绘制各含水层初始流场、末刻流场的基础上，模拟计算出剖分图上各含水层全部节点的初始水位值、末刻水位值，作为数值法水量计算程序所用的初始流场、及对照末刻计算流场的实测流场。

7.模型识别

把所确定的模型识别起始日期到结束日期之间的源汇项数据按时段代入计算程序、把各区的参数初值作为模型的参数初值代入计算程序，运行程序，通过不断调试模型参数得到的拟合孔的计算历时曲线与实测的历时曲线进行拟合，调试出一组仿真度最高的模型参数。

8.模型检验

把所确定的模型检验起始日期到结束日期之间的源汇项数据按时段代入计算程序、把模型识别调试出的模型参数代入计算程序，主要看模型识别结束日期之后拟合孔的计算历时曲线，与实测的历时曲线的吻合程度，来检验模型识别阶段调试出的模型参数是否经得起“外推”的检验；如果“外推”段的拟合曲线误差太大，则应推倒识别阶段的模型参数，重新进行模型识别、模型检验，直到满足误差要求为止。

9.模型预报

经过模型识别、通过模型检验后确定了模型参数，即最终确定了数值模型，即可用其来进行模型预报。进行模型预报时，一般根据供水需要设计多个可供选择的开采方案，分别代入数值模型，来预测未来数十年（城市供水水源地一般不少于30a）后的流场，然后，通过对比、分析，选择其中较合理的开采方案作为水源地建设、施工设计的水文地质依据。

本章以下几节分别介绍不同水文地质条件的几种二维流有限元水量计算问题。

一、太原盆地水文地质参数计算

水文地质参数的选取直接影响着地下水资源计算量的大小和可信度，研究水文地质参数具有十分重要的意义。本次相关的水文地质参数主要有降水入渗补给地下水系数（α）、潜水蒸发极限深度（L）、蒸发强度（ε）、灌溉回渗地下水系数（β）、疏干给水度（μ）、导水系数（T）、d性储水系数（s）、渗透系数（K）、河流渗漏补给系数、渠系渗漏补给系数等。

（一）降水入渗补给地下水系数（α）

影响降水对地下水的补给量的因素很多，主要有地形、包气带岩性及结构、地下水位埋深、降水特征及土壤前期含水量等。

降水入渗补给系数为降水入渗补给地下水量与降水量之比值。年降水入渗补给系数为年内所有场次降水对地下水入渗补给量总和与年降水总量的比值，其表达式为:

山西六大盆地地下水资源及其环境问题调查评价

式中:α年是年降水入渗补给系数；pri是场次降水入渗补给量，mm；P是年降水量，mm；n是年降水场次数。

用长期动态观测孔求取年降水入渗系数的计算方法:

山西六大盆地地下水资源及其环境问题调查评价

式中:μ∑Δh次是年内各次降水入渗补给地下水量之和；P年是年降水量；Δh次是某次降水引起的地下水位升幅值。

根据动态资料分析计算，在前人试验的基础上，综合考虑各方面的因素，给出盆地区降水入渗补给地下水系数（详见第四章）。

（二）地下水蒸发极限深度（L）、蒸发强度（ε）

蒸发极限深度就是指浅层水停止蒸发或蒸发量相当微弱时，浅层水位埋深值。蒸发强度就是在极限蒸发深度以上，单位时间浅层水的蒸发量。

影响地下水蒸发的主要因素是地下水位埋深、包气带岩性和水面蒸发强度等。

理论上，当水位埋深处于蒸发极限深度时，地下水在无补给、无开采的条件下，动态曲线近于平直。

地下水蒸发极限深度（L）

蒸发极限深度通常采用迭代法、试算法和经验公式计算（L），公式如下:

迭代法:

试算法:

经验公式法:

式中:ΔT1、ΔT2为计算时段，d；H1、H2、H3为时段内水位埋深，m；Z1、Z2为时段内水面蒸发强度，m/d；

经计算，太原盆地孔隙水区不同岩性的蒸发极限深度依包气带岩性不同分别为:亚砂、亚粘土互层为3.5m，亚砂土为4.0m，粉细砂、亚砂土互层为4.5m。

地下水蒸发强度

计算公式:

式中:Z0是液面蒸发强度，mm/d；ΔH是浅层水降落间段的平均水位埋深，mm；Z是蒸发强度，mm/d。

由本区浅层水水位埋深图（详见第四章）可看出，水位埋深小于4m的区域在北部太原市和南部平遥、介休一带，根据上式计算太原、平遥、介休等地的地下水蒸发强度见表3-1。

表3-1 太原盆地孔隙水区地下水蒸发强度

（三）灌溉回渗地下水系数（β）

是指田间灌溉补给地下水的量与灌溉总量的比值。影响灌溉回渗系数和因素主要有岩性、水位埋深、土壤含水率、灌溉定额等多种。

计算公式:

式中:μ是给水度；Δh是由灌溉引起的地下水位平均升高值，m；Q是灌溉水量，m3；F是面积，m2。

本次工作在盆地太原市小店区郜村、汾阳市贾家庄镇东马寨村和榆次市杨盘等3个地方布置了3组灌溉入渗试验，地表岩性郜村为粉质粘土、东马寨上部为粉质粘土，下部为粉土，杨盘为粉土，化验室给水度试验结果分别为0.195、0.11、0.143。郜村在37m×37m的面积上布置10眼观测孔，水位埋深1.2～1.3m，累计灌溉水量160m3，10个孔平均水位上升值为0.1912m，根据上式计算得灌溉入渗地下水系数为0.32；东马寨村水位埋深1.95～2.44m，在26m×26m的面积上布置10眼观测孔，灌溉水量60m3，观测孔平均水位上升值为0.465m，计算得灌溉入渗地下水系数为0.58；杨盘布3个观测孔，水位埋深5.76～6.01m，灌溉面积100m2，灌溉水量100m3，平均水位上升高度为0.27m，计算得灌溉入渗系数为0.039。

从以上试验数据可以看出，不同水位埋深、不同岩性地区灌溉入渗系数有很大区别。综合考虑各种因素，灌溉回渗地下水系数选用值见表3-2。

表3-2 灌溉回渗地下水系数

（四）d性贮水系数S、导水系数T、给水度μ、渗透系数K

盆地区大部分地区都进行过1∶5万比例尺的农田供水水文地质勘查，做过大量单孔和多孔抽水试验，本次在文水文倚、汾阳等5地分别作了5组抽水试验，用非稳定流公式，降深-时间半对数法计算结果如下:文倚导水系数T＝1983.59～2181.95m2/d，渗透系数K＝32.19～35.4m/d，d性贮水系数S＝1.79×10-3；汾阳县贾家庄镇东马寨村抽水试验求得导水系数T＝325.84～376.5m2/d，渗透系数K＝5.65～6.53m/d。结合以往本区的工作成果，给出太原盆地浅层孔隙潜水和中深层孔隙承压水水文地质参数，详见参数分区图3-13和参数分区表3-3。

表3-3 太原盆地中深层孔隙承压水及浅层孔隙潜水参数分区

图3-13 太原盆地参数计算分区图

二、大同盆地水文地质参数计算

由本区浅层水2004年水位埋深图可看出，水位埋深小于4m的区域主要分布于盆地中部冲积平原区，盆地南部怀仁、山阴、应县、朔州分布面积较大。根据计算和以往试验资料，本区蒸发强度确定值见下表（表3-4）。

表3-4 大同盆地孔隙水区地下水蒸发强度

据“山西省雁同小经济区水资源评价、供需平衡研究报告”中搜集的本区灌溉回渗试验数据取得不同水位埋深、不同岩性、不同灌溉定额的灌溉回渗系数，灌溉回渗系数选定值见表3-5。

盆地区大部分地区都进行过1/5万比例尺的农田供水水文地质勘查，做过大量单孔和多孔抽水试验。本次工作搜集本区以往抽水试验孔117个，本次在大同县党留庄乡、怀仁县金沙滩镇、怀仁县新发村、怀仁县榆林村、山阴县张庄乡、朔州市城区沙塄乡等6地分别作了6组抽水试验，采用AquiferTest计算程序，非稳定流方法计算，本次抽水孔具体情况和计算结果见表3-6和表3-7 。

表3-5 灌溉回渗地下水系数

表3-6 大同盆地本次抽水试验数据统计

表3-7 大同盆地本次抽水试验计算成果表

结合以往本区的工作成果，给出大同盆地浅层孔隙潜水和中深层孔隙承压水水文地质参数，详见参数分区图3-14、图3-15和参数分区表3-8、表3-9 。

图3-14 大同盆地降水入渗系数分区图

图3-15 大同盆地浅层、中深层孔隙水参数分区图

表3-8 大同盆地浅层孔隙潜水参数分区表

续表

表3-9 大同盆地中深层孔隙承压水参数分区

三、忻州盆地

忻州盆地地下水资源较为丰富，开采条件优越，20世纪70年代之前地下水开采规模较小；70年代初至80年代末随着农业灌溉的普及，工业生产的发展和城市规模的扩大，地下水开采量迅速增加。开采对象以浅层水为主，造成浅层水水位普遍有所下降（但下降幅度不大）。从20世纪90年代至今，虽然地下水开采量具有逐年增大的趋势，但增加幅度较小，且中层井数量逐渐增多，形成了浅层水、中层水混合开采的新模式，地下水位总体处于动态平衡状态。受地下水人工开采的影响，降水入渗系数及导水系数等水文地质参数发生了一定程度的变化。

区内降水入渗系数的变化除了与年降水量及降水特征有关外，主要与浅层地下水位埋深关系较为密切。已有资料表明，在山前倾斜平原区，浅层水位埋深一般大于7m，因水位下降使降水入渗系数发生了不同程度的减小。在冲积平原区浅层水位埋深一般小于7m，水位下降的结果引起了降水入渗系数有所增大。不同地貌单元降水入渗系数的变化见第五章。

从20世纪70年代以来，区内含水层的导水系数发生了较为明显的减小，主要体现在因浅层地下水位下降，使浅层含水层上部处于疏干状态，含水层厚度减小，直接导到导水系数减小。因浅层水水位下降幅度不同，导水系数减小的程度也存在差异，从本次地下水侧向补给量计算断面附近的井孔资料分析，含水层厚度一般减小了3～6m，导水系数由70年代中期的60～250m2/d，减少到目前的50～200m2/d左右。

忻州盆地给水度根据不同地貌单元含水层岩性、分选性及富水性综合确定见表3-10及图3-16 。

表3-10 忻州盆地浅层含水层给水度分区

图3-16 忻州盆地给水度分区图

四、临汾盆地

经过搜集以往资料，调查和计算确定临汾盆地降水入渗系数见表3-11。临汾盆地渗透系数及给水度分区见图3-17，表3-12。

表3-11 临汾盆地平原区降水入渗系数统计

图3-17 研究区渗透系数及给水度分区图

表3-12 临汾盆地参数分区表

五、运城盆地

运城盆地地下水长观网建站年代较远，积累了大量的地下水位监测资料，且经过多次的地质、水文地质勘察、地下水资源评价工作，取得了大量的降水入渗值，参考前人综合成果，结合目前包气带岩性、地下水位埋深，给出运城盆地降水入渗补给系数，见表3-13。

表3-13 运城盆地平原区降水入渗系数统计

渠系有效利用系数除受岩性、地下水埋深影响外，还与渠道衬砌程度有关。修正系数r为实际入渗补给地下水量与渠系损失水量Q损的比值，是反映渠道在输水过程中消耗于湿润土壤和侵润带蒸散损失量的一个参数，它受渠道输水时间、渠床土质及有无衬砌、地下水埋深等因素的影响。一般通过渠道放水试验获得。本次评价主要参考运城市水利局相关试验成果，见表3-14。

表3-14 运城盆地万亩以上灌区η、r、m值统计

灌溉回归补给系数β值与岩性、植被、地下水埋深及灌溉定额有关，一般通过灌溉入渗试验求得，本次评价主要参照运城市水利部门资料综合确定，详见表3-15。

表3-15 运城盆地灌溉回归系数β取值

河道渗漏补给系数是河道渗漏补给地下水量与河道来水量的比值。其值大小与河床下垫面岩性、流量、地下水位埋深及渗漏段长度有关。运城盆地沿中条山前发育数条季节性河流，河床下垫面主要为砂卵砾石，当洪雨季节，地表河床水位远高于地下水位，为地表水的入渗造就了十分便利的条件。根据河道渗漏资料，可建立如下数学模型:

山西六大盆地地下水资源及其环境问题调查评价

式中:m河是河道渗漏补给系数；A是计算系数，A＝（1-λ）×（1-φ）L，φ是单位千米损失率；L是河道渗漏长，km，Q径是河道来水量，m3/s。

据运城市水利部门研究成果，A值约为0.090。

含水层的渗透系数主要由野外抽水试验通过稳定流及非稳定流计算公式求得，各勘探部门在运城盆地先后进行过各种勘察，进行了大量的抽水试验工作，积累了丰富的资料，参考本次抽水试验成果对以往参数进行了修正，取值结果见表3-16 。

表3-16 运城盆地松散岩类K值选定表

降雨入渗补给系数在同岩性、同降雨量情况下，随地下水位埋深的增大，降雨入渗补给系数会达到一个最大值之后趋于减少或变为常数。运城盆地北部的峨嵋台塬及闻喜北塬，其地下水位埋藏深，地表主要以黄土类为主，降水入渗主要依靠黄土垂直节理裂隙及“流海缝”以“活塞式”注入地下，多年来其降水入渗系数基本为常量，经用动态分析法计算其降水入渗系数在0.108～0.11间；在盆地中部的冲湖积平原区，其地表岩性主要以Qp3＋Qh冲湖积相的亚砂土、亚粘土、粉细砂为主，由于开采强烈，区域水位严重下降，地表数米至几十米内均为饱气带，为降水入渗准备了调蓄空间，加强了降水向地下水的转化。根据盆地地下水长观孔资料及次降雨资料，计算出盆地冲湖积平原地带，降水入渗系数在0.1～0.162之间，总体上上游大于下游。而在东部及南部的山前倾斜平原区，地下水位埋深一般大于5m、乃至几十米，地表岩性大多为亚砂土及亚粘土，尤其是在一些沟口附近，从地表往下几十米范围内为干砂卵砾石，一般降雨基本上不产生地表径流，这无疑加大了降水的转化。据相关资料计算，降水入渗系数高达0.21～0.30。因过去所做的工作不系统，没有对降雨入渗系数进行系统分类，不便比较，但根据运城盆地饱气带岩性、地下水变动情况，除峨嵋台塬及黄土丘陵区变化不大外，其他地区降雨入渗系数无疑有增大趋势。

盆地内抽水井的含水层，大多为数个含水层混合开采。现根据本次抽水计算值，对历次研究成果中的K值加以修正，得出运城盆地各个地貌单元的渗透系数。总体来说，黄河岸边低阶地区K值最大为11.3～14.6m/d，中条山山前倾斜平原次之，为5.45～6.12m/d，最次为闻喜北垣K＝1.10m/d左右。

根据地貌单元、含水层岩性、地下水水力特征及各参数特征，将运城盆地划分为10个参数分区，见表3-17及图3-18。

表3-17 运城盆地水文地质参数分区

六、长治盆地

根据水文地质条件，长治盆地参数分区见图3-19，表3-18 。

图3-18 运城盆地水文地质参数分区表

图3-19 长治盆地参数分区图

表3-18 长治盆地浅层孔隙潜水参数分区

（一）降水入渗补给系数变化

根据《太原市地下水资源评价报告》研究成果，盆地区亚砂土、极细砂、细砂的降水入渗系数随着地下水位埋深的增大而增大，当水位埋深超过一定值以后，降水入渗系数开始趋于稳定；降水量越大，降水入渗系数在相同的岩性和地下水位埋深条件下也越大。对于亚砂土、极细砂、细砂在相同水位埋深和降水情况下，细砂的降水入渗系数>极细砂的>亚砂土的。总体来说，颗粒越粗，降水入渗系数也越大。

α随降水量的变化，非饱和带在降水入渗补给地下水过程中起调节作用，降水入渗补给过程要滞后于降水过程，其滞后时间的长短、特征与非饱和带的重力水蓄水库容关系密切，地下水埋深越大，其蓄水库容也越大，调节能力也越强，滞后现象也越明显。

在亚砂土、极细砂和细砂3种岩性中，降水量相等时，降水入渗系数从大到小的顺序为细砂、极细砂、亚砂土。场次降水量的影响表现为α次先是随着降水量的增大而变大，当降水量超过一定数值后，α次反而呈减少趋势，这个降水量即是最佳降水量。α年与α次有相同的规律性，从入渗机制分析，α年也存在最佳年降水量。

当地下水埋深为零时，降水入渗补给系数亦为零，然后随埋深的增加由小变大；当地下水埋深到达某一定值时，降水入渗补给系数达到最大值即最佳降水入渗补给系数，并由此随埋深的增加由大到小，到达一定的埋深时，趋于定值。地下水埋深对降水入渗补给系数的影响，可从3方面来说明。

埋深反映了蓄水库容的大小。当埋深为零时，即蓄水库容为零，这时无论降水量多大，均无入渗补给的可能。当埋深增加时，地下水库得到了降水入渗补给量，此时降水入渗补给系数大于零，降水入渗补给系数随埋深的增加而增大。当地下水达到最佳埋深时，其对应的降水入渗补给系数为最佳降水入渗补给系数，原因是由于条件一致的地区中的依次降水，其入渗补给量随地下水埋深的变化必存在一个最大值。当地下水埋深较小时，由于地下水蓄水库容较小，形成蓄满产流，不能使降水全部入渗；当地下水埋深再增大时，则损失较最佳埋深为大，故降水入渗补给系数随埋深的增加而减小。对于不同级别的降水量，α最大值出现的地下水位埋深区域也不同。最佳埋深与岩性和降水量有关。

地下水埋深在某种程度上反映了土壤水分的多少。土壤水垂直分布大体可概化为3种状况。第1种情况是地下水埋深较小，毛管上升水总能到达地表；第2种情况是地下水埋深较大时，毛管上升水无法到达地表；第3种情况是地下水埋深介于两者之间，在此埋深内，由于地下水位是升降变化，毛管上升水有时达到地表，有时达不到地表。这3种情况将对降水入渗补给量有不同的影响。第1种情况，降水一开始，水即可通过毛管在重力作用下迅速向下移动，地下水位在降水开始后很快上升。第2种情况，降水首先应满足土壤缺水的需要，而后在重力作用下通过空隙下渗补给地下水。其渗漏途径较第1种情况长，入渗方式也有差异。

图3-20 渗透系数与深度关系图

不同地下水位埋深条件对降水入渗补给系数取值的影响。盆地太谷均衡实验场的水分势能实验最大深度为8.2m，有观测点41个。多年资料的分析结果表明，土壤水分势能变化从地面往下可分为3个变化带———剧烈变化带、交替变化带和稳定带，剧烈变化带埋深为0～1.1m，土壤水分势能变幅大于200×133Pa；交替变化带埋深1.1～3.6m，土壤水分势能变幅大于（100～200）×133Pa之间；埋深3.6m以下为稳定带，其土壤水分势能变幅小于100×133Pa，其中埋深在4.5～5.0m以下的稳定特性更为明显，其土壤水分势能的变幅一般不超过50×133Pa，其土壤水分全年为下渗状态。表明埋深在5.0m以下为稳定入渗补给，反映在降水入渗补给系数上随埋深增加，α年将趋于稳定，故当埋深大于5.0m时，α年值可取定值，不再随埋深而变化。原因是地下水埋深已到达或超过地下水极限埋深，损失趋于定值，水分不向上运动，必然向下运动，故形成了降水入渗补给系数随地下水埋深变化的稳定值。

（二）渗透系数变化

孔隙含水介质的渗透能力不仅取决于粒径大小、颗粒级配、胶结程度，还与其埋深有关。同一岩性的孔隙含水介质，随着深度的增加，介质被压密，渗透系数会减小。

根据河北平原山前冲洪积扇扇顶区数百个钻孔资料的统计，各种含水介质的渗透系数随埋深增加呈指数衰减，部分深层不同岩性渗透系数随埋深的变化规律参考下述经验公式:

岩性为卵砾石时，渗透系数与埋深关系式:

K＝K0e-0.0131h R＝0.877

岩性为砂砾石时，渗透系数与埋深关系式:

K＝K0e-0.0116h R＝0.869

岩性为中粗砂时，渗透系数与埋深关系式:

K＝K0e-0.0057h R＝0.896

K为埋深处的渗透系数；K0为地表浅层的渗透系数；h为埋深；R为相关系数。

因此，对于同一种岩性，其渗透系数大小与深度有关（图3-20）。

数值模型模拟计算方法适用于非均质性、各向异性的复杂地下水系统，包括存在越流和具有不规则形状各类边界条件等情况。但是该方法对资料的要求比较严格，要求研究程度较高和资料较丰富。应用数值模型方法的一般程式为：①气象、水文、水文地质资料分析→②水文地质概念模型概化→数学模型建立（水动力方程和定解条件）→③选择计算程序→④模型设计→⑤模型识别和检验→⑥计算模拟。

一、水文地质概念模型建立

在对黑河流域地下水系统做了全面、深入分析的基础上，根据研究目的，对地下水系统的组成要素和相互关系作出合理的简化和假设，并且用文字、框图、平面图、剖面图等形式把系统再现出来，即为地下水系统概念模型。

（一）地下水系统空间结构概化与边界确定

1.图式表示地下水系统空间结构

根据黑河流域水文地质图和水文地质剖面图，梳理和划分主要含水层、隔水层与弱透水层，阐明它们的产状、分布范围和厚度等，确定透水、阻水等断层属性。分析地下水系统的各类等值线图，包括第四系基底埋深等值线图、地下水水头等值线图、含水层顶底板高程等值线图、含水层和隔水层的厚度等值线图等。

2.确定地下水系统边界

地下水系统的边界，包括自然边界（固定边界）和水力边界（可移动边界）。自然边界包括不透水岩层、不透水断层或断裂带、较大的地表水体等；水力边界包括地下水分水岭和地下水流线等。

数值模型模拟研究，其对象的底界一般为不透水岩层。侧向边界可以是自然边界，也可以是水力边界或无穷远边界（边界水头或流量不受输入条件的影响）。模拟顶界，对于承压水系统而言，一般为不透水边界或越流边界，对于潜水系统一般采用大气边界（蒸发和入渗）。地下水系统内部边界包括零流量边界（不透水岩体）和流量边界（河流、湖泊或水库的渗流带）等。

3.水文地质参数

水文地质参数是数值模型模拟研究的灵魂，一般包括含水层组的渗透系数、导水系数、给水度、储水率、储水系数、孔隙度、垂向渗透系数和越流系数，以及包气带的降水入渗系数、河道渗漏系数、井灌回归系数、田间与渠道渗漏系数、潜水蒸发系数和陆面蒸发系数等。

确定降水入渗补给系数、灌溉渗漏系数、蒸发系数等方法，有水文分析法（降水量、河流径流量曲线、地下水水头动态曲线等）、直接试验法（地渗仪、张力计、同位素示踪等）、计算法（氯质量平衡法、非饱和模型法等）、经验公式法和ZFP零通量面实测法等。

（二）地下水流系统概化

对地下水流系统进行概化，包括确定地下水的基本流向、地下水补给要素组成、排泄模式、地下水与地表水之间转化关系、不同层位含水层之间水力关系等。主要依据有，地下水水头等值线图、水化学信息、同位素信息、地下水温度信息和水位动态曲线等。

根据地下水流状态及其特征，确定所研究的地下水流系统具体属性，例如稳定流或非稳定流，一维流、二维流、准三维流或三维流等。

（三）模型输入量计算

降水入渗、地表水入渗（河渠）、地下水侧向流入、灌溉入渗、蒸发蒸腾、泉水排泄、基流排泄、地下水侧向流出、开采等。

二、建立数学模型

根据建立的水文地质概念模型，选择适宜数学模型。一般由描述地下水运动规律的偏微分方程和反映地下水系统边界条件及初始条件的定解条件组成。

非均质承压水三维非稳定流偏微分方程为

西北内陆黑河流域水循环与地下水形成演化模式

非均质无压水三维非稳定流偏微分方程有下列几种情况：第一类边界条件（狄利克雷边界）为

西北内陆黑河流域水循环与地下水形成演化模式

第二类边界条件（纽曼边界）为

西北内陆黑河流域水循环与地下水形成演化模式

初始条件为

西北内陆黑河流域水循环与地下水形成演化模式

三、计算程序、模型设计与识别

（一）计算程序与模型设计

计算程序分为一维流、二维流、准三维流或三维流模型，以及对均质、非均质、各向同性或各向异性和对不同输入项的处理能力。目前，可供软件有MODFLOW、FEWFLOW、PM、GMS、GWVISTA、MODME、PM等，它们多为有限差分法和有限元法。模型设计，包括网格剖分（规则剖分或不规则剖分、三角剖分或矩形剖分）、选择时间步长（试算法）、设置模型边界、设置初始条件、数据输入（降水入渗速率、田间灌溉入渗速率、蒸发速率、水井位置及开采或回灌强度、地下水与地表水相互作用的时空分布、泉的时空分布、边界水位或边界流量、观测井位置及观测水位等）。

（二）模型识别与检验

1.模型识别

模型识别亦称反演问题，即利用实测地下水动态资料和抽水试验资料，反求水文地质参数或源汇项和定解条件的过程。模型识别是为了解决选用的偏微分方程是否合适问题，确定模型中的水文地质参数和源汇项及定解条件，从而建立一个能再现地下水系统实际功能（水头或浓度）的模拟模型。模型识别一般采用试估-校正法。就是选择一合适的时段，根据水文地质条件和经验数据估算一组水文地质参数输入模型，利用所选时段的输入输出数据，求解模型。然后以模型计算结果与实测结果比较，如果拟和结果不符合精度要求，适当调整参数，重复上述过程，直到符合精度要求为止。也可以采用试估-校正法与最优化方法相结合的方法。首先用试估-校正法粗调，然后用最优化方法细调，即用最优化方法求得一组最佳的参数值，使得计算水头值与观测值之间的差值在给定的约束条件下，达到极小。

模型识别的结果具有多解性。要识别的参数数目应少于总数据数目。也就是说必须要有已知量。已知量愈多，反求的参数愈精确，由此建立的模型的适用性就愈好。正因为模型识别结果的多解性，所以对于同一个问题，不同的人所求得的参数组合不同，甚至同一个人在不同的时间所求得的参数也不同。显然，模型识别的参数不一定是含水层所固有的参数。因此，有人称模型识别的参数为“模型参数”，以示区别。尽管模型参数不能完全反映实际系统的参数，但是模型参数有其特殊作用，它能够使得数学模型在行为和功能上代替实际的地下水系统，成为地下水系统的“复制品”。

2.模型检验

为了检验识别后的模型的可靠性，需要采用同一系统的另一时段的数据资料输入模型进行检验。如果计算结果符合实际资料，则可以说明模型能真实反映实际系统。需要指出的是，在模型识别和模型检验阶段所用的两组数据资料，必须是相对独立的不同时间段的资料。

模型灵敏度分析的目的是了解参数变化对计算结果的影响，同时识别重要参数。灵敏度分析一般在模型识别之前进行，也可以在模型识别之后进行。

选取要分析的一个参数（θ），然后固定其余参数，改变θ的数值分析计算结果。这时计算水头（g）就是θ的函数，即g=f（θ）。则有如下定义：在θ=θ0附近，水头变量g（θ）相对于原值g*（θ）的变化率和参数θ相对于θ0的变化率之比称为水头对参数θ的灵敏度，以下式表示：

西北内陆黑河流域水循环与地下水形成演化模式

四、黑河流域模拟区水文地质条件概化

地下水数量转化研究的数值模型模拟区，选择了张掖盆地和酒泉东盆地，包括张掖、临泽、高台的所有灌区和民乐及山丹的个别灌区，还有肃南县明花区，面积近9000 km2。

数值模拟区是只有侧向流入而没有侧向流出的山间断陷盆地，其间充填了巨厚的松散沉积物，构成赋存地下水的天然场所，为连续和统一的第四纪含水岩系综合体，周边山体为天然的地质边界。在张掖盆地，地下水自南东向北西运动，泄于黑河干流而流出区外。西部酒泉东盆地，地下水由南西向北东运动，榆木山至高台县城一线为两盆地天然汇水线。

数值模拟区地下水的主要补给来源是河水（含雨洪水）、渠系引水和田间灌溉水的垂直入渗，而泉水溢出、蒸发和人工开采是主要排泄方式。

据均衡计算结果，1999年区内补给量为11.94×108 m3，排泄量为14.09×108 m3，均衡差为-2.15×108 m3，数值模拟区处于负均衡状态，地下水水位呈下降态势。

数值模拟区周边皆为二类流量边界。山区边界沿山前大断裂分布，流入量主要为基岩裂隙水侧向流入和沟谷潜流。东部民乐、山丹断面和西部明花区断面为区外侧向流入量，利用断面法求得。南部新坝-红崖子隐伏断层使地下水流不连续，作为该段边界，概化的水文地质模型如图5-1。

五、数学模型概化

数值模拟区南半部为潜水、北半部为承压水，适宜采用潜水-承压水数学模型。但是各灌区开采地下水的程度不同，一些地带已将潜水与承压水连通，承压水头与潜水水位动态变化具有一致性。因此，将模型概化为非均质各向同性二维流潜水模型。鉴于区域面积大，地下水水位年变幅小，与含水层厚度相比可忽略，所以用导水系数（T）近似代替渗透系数（K）与含水层厚度（H）之积。

数学模型及定解条件如下：

图5-1 黑河流域模拟区水文地质模型概化图

西北内陆黑河流域水循环与地下水形成演化模式

式中：T——含水层导水系数（m2/d）；

μ——含水层给水度（无量纲）；

Wb——各项补给项强度之和（m3/km2·d）；

Wp——各项排泄项强度之和（m3/km2·d）；

q——流量边界单宽流量（m3/km2·d）；

Γ2——流量边界代号；

n——边界上的内法线方向。

采用线性插值，伽辽金有限元法解上述方程组，见程序框图（图5-2）。

图5-2 黑河流域数值模型模拟程序求解流程

六、定 解 条 件

（一）初始条件

以1999年水位统测结果为基础，结合地下水动态长观资料，绘制1月份等水位线图为初始流场。采用三角剖分法将计算区剖分成1421个单元，799个结点。其中内结点624个，边界点175个。水位观测点33个，均分布于结点上（图5-3）。同时尽量把结点布置在概化的灌区边界上。

（二）计算时段

以1999年元月初至12月末每个自然月实际天数为时段长度，全年共分12个时段。

（三）水文地质参数

根据黑河勘察报告研究成果，数值模拟区参数取值范围T值为100～6500 m2/d，μ值为0.1～0.25之间。参数分区以灌区为基础，按不同埋深划分。

（四）源汇项

计算区地下水主要靠河水、渠系引水、灌溉水、降水凝结水入渗及边界流入补给。消耗于蒸发蒸腾、泉水溢出和人工开采。有关参数的选取，主要依据黑河报告和各县水利部门研究成果，补给量与排泄量通过水量均衡方法计算求得。

由于数值模拟区范围较大，而且区内农业发达、干支渠密布，沿主要河流（黑河）引水口众多，所能收集到的水文和水利资料有限，所以剖分不宜过细，可将河水（含雨洪）、渠系水、灌溉水和降凝水入渗及人工开采处理为面状量，把各灌区不同埋深均衡计算结果以单位面状量进入模型，补给项为正，排泄项为负。非灌溉期（1～3月，10～12月）的渠系水和灌溉水入渗及人工开采量强度为0，灌溉期（4～9月）摊分全年入渗量。

图5-3 黑河流域数值计算区剖分图

1999年河水入渗量占当年黑河（莺落峡）径流量的32%，每月径流量占全年径流量的比例分配到12个时段。降水、蒸发强度按各月份所占全年比值分配到12个时段。1～3月和10～12月的降水为0，4～6月降水占30%，7～9月降水占70%。按地下水水位不同埋深，计算蒸发量，其中1～3月占13%，4～6月占41%，7～9月占35%，10～12月占11%。

泉水溢出带均分布于细土平原、地下水水位埋深小于3.5m的地带，各泉沟及黑河河床地下水水位高于河床标高，实际为线状量。但是因剖分单元较大，无法准确描述，所以将线状量处理成面状量，假设地下水水位埋深小于3.5m带为泉水溢出带，具体做法将所有结点地面高程减去3.5m，于是该区地下水水位埋深值为负。将1999年泉水溢出量除以该区面积，再除以平均水头差1.5m，获得单位水头差条件下泉水溢出强度，引入模型。然后根据各时段水头变化，获得不同时段的泉水溢出量。

数值模拟区边界为透水边界或弱透边界，均给出单宽流量，全年一致，不再按时段划分。

七、数 模 结 果

按上述补给与排泄要素及其参数，采用观测点的地下水水位拟合，对1999年实施模型进行识别。

区内共有观测点33个，集中在张掖、临泽、高台的细土平原带。在调参过程中，不断缩小拟合点误差，兼顾初始流场与计算流场形态一致，并且每个节点水位偏差不宜过大。调参结果，数值模拟区共有60个参数分区，如图5-4和表5-2所示。观测点拟合结果如图5-5和图5-6所示，地下水流场拟合情况如图5-7所示。

图5-4 黑河流域数值模拟参数分区图

表5-2 黑河流域模型采用的有关水文地质参数

---