数值模型模拟计算方法适用于非均质性、各向异性的复杂地下水系统，包括存在越流和具有不规则形状各类边界条件等情况。但是该方法对资料的要求比较严格，要求研究程度较高和资料较丰富。应用数值模型方法的一般程式为：①气象、水文、水文地质资料分析→②水文地质概念模型概化→数学模型建立（水动力方程和定解条件）→③选择计算程序→④模型设计→⑤模型识别和检验→⑥计算模拟。
一、水文地质概念模型建立
在对黑河流域地下水系统做了全面、深入分析的基础上，根据研究目的，对地下水系统的组成要素和相互关系作出合理的简化和假设，并且用文字、框图、平面图、剖面图等形式把系统再现出来，即为地下水系统概念模型。
（一）地下水系统空间结构概化与边界确定
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 黑河流域模型采用的有关水文地质参数
渠系渗漏修正公式
一、均衡区的划分及时段选择
依据泾惠渠灌区地质构造、地貌单元进行水文地质分区，将灌区分为4个区，7个亚区。泾、渭河漫滩及一级阶地的第四系全新统含水岩组区；泾河二级阶地第四系全新统含水岩组区；泾河三级阶地、渭河二级阶地第四系上更新统含水岩组区；黄士台源第四系中更新统含水岩组区。泾河漫滩及泾河一级阶地强富水亚区；渭河漫滩及渭河一级阶地强富水亚区；泾河二级阶地清峪河以南富水亚区；泾河二级阶地清峪河以北富水亚区；泾河三级阶地中等富水亚区；渭河二级阶地中等富水亚区；黄士台源弱富水亚区（陕西省泾惠渠灌区地下水调查研究组，1983）。
均衡时段的确定，根据灌区灌溉年度即11月1日起至次年10月底为均衡计算时段。各灌季的时段分为：冬灌11月1日至2月底；春灌3月1日至5月底；夏灌6月1日至8月底。
二、均衡方程的建立
依据水均衡原理，结合灌区潜水的补给、径流、排泄条件，建立如下的水均衡方程（地矿部水文地质工程地质技术方法研究所，1978）：
灌区农业节水对地下水空间分布影响及模拟
其中：
Q补＝Q田补+Q雨补+Q渠补+Q侧补+Q开补
Q排＝Q开采+E+Q侧排+Q用
式中：Q补为地下水总补给量，m3/d；Q排为地下水总排泄量，m3/d；μ为水位变动带给水度；F为均衡区面积，km2；?t为均衡时间段长，d；?H为与 ?t对应的水位变幅，m；Q田补为渠灌田间入渗补给量；Q雨补为降水入渗补给量；Q渠补为渠系渗漏补给量；Q侧补为侧向径流补给量；Q开补为井灌开采回归补给量；Q开采为地下水开采量；E为潜水蒸发量；Q侧排为侧向排泄量；Q用为人畜工业用水开采量。
三、补给量计算
1.降水入渗补给量（Q雨补）
降水入渗量的多少，主要受地层岩性等地质条件影响，根据灌区气象站2010年降水量408.0mm及各典型年的降水量，将其代入下式中，计算出大气降水入渗补给量。
Q雨补＝∑ai·Pi·Fi （7-8）
式中：ai为不同区段的降水入渗系数值；pi为不同频率下的降水总量（mm）；Fi为不同区段的面积，km2。灌区降水入渗补给量见表7-5。
表7-5 灌区降水入渗补给量 Table7-5 precipitation infiltration recharge in Jinghui Canal Irrigation District
2.地下水侧向流入量（Q侧补）
根据所选的断面位置、断面长度、含水层平均厚度、平均水力坡度、平均渗透系数，利用下式计算地下水侧向径流补给量。
Q侧补＝K·I·B·H （7-9）
式中：K为含水层的渗透系数，m/d；I为断面处的水力坡度；B为断面宽度，km；H为潜水含水层厚度，m。2001年侧向补给量见表7-6。
表7-6 灌区2005年侧向补给量 Table7-6 Lateral recharge in Jinghui Canal Irrigation District in 2005
3.渠灌田间入渗补给量（Q田补）
本次将斗、农、毛三级渠道的渗漏补给量计入渠灌田间入渗补给量。渠灌田间入渗补给量可利用下式计算：
Q渠灌＝β渠·Q渠田 （7-10）
式中：Q渠灌为渠灌田间入渗补给量，104m3；β渠为渠灌田间入渗补给系数（无因次）；Q渠田为渠灌水进入田间的水量，104m3（应用斗渠渠首引水量）。
利用式（7-10）计算多年平均渠灌田间入渗补给量见表7-7，Q渠田采用1990～2010年期间的多年平均值，β渠采用近期地下水埋深和灌溉定额条件下的分析成果。
表7-7 灌区渠灌田间入渗补给量 Table7-7 Canal irrigation field infiltration recharge in Jinghui Canal Irrigation District
4.井灌回归补给量（Q开补）
根据灌区实际统计井灌面积计算井灌回归补给量，井灌回归补给系数β井统一取0.17。
5.渠系渗漏补给量（Q渠补）
本次只计算干、支两级渠道的渗漏补给量，利用干支渠计算时段引水量和渠系渗漏损失系数计算，计算公式：
Q渠补＝m·Q渠首引 （7-11）
式中，Q渠首引为渠首引水量，104m3；m为渠系渗漏补给系数（无因次）。利用渠系渗漏补给系数法，即利用式（7-10）计算多年平均渠系渗漏补给量Q渠补时，相关计算参数应采用2001～2009年期间的渠系渗漏补给量（表7-8）。
表7-8 灌区渠系渗漏补给量 Table7-8 Canal system leakage recharge in Jinghui Canal Irrigation District
四、排泄量计算
1.地下水开采量（Q开采）
地下水开采量指灌溉开采的地下水量，根据灌区实际调查统计的井灌面积与灌水定额来计算。
2.潜水蒸发量（E）
潜水蒸发经验公式用修正后的阿维里扬诺夫（C.φ.AвepъянОв）公式：
灌区农业节水对地下水空间分布影响及模拟
式中：Z0为极限埋深，m，即潜水停止蒸发时的地下水埋深；n为经验指数，一般为1.0～2.0，应通过分析，合理选用；k为作物修正系数；Z为潜水埋深，m；E、E0分别为潜水蒸发量和水面蒸发量，mm。由于本区地下水位的埋深均＞7m，因此忽略不计潜水面的蒸发量。
3.侧向排泄量（Q侧排）
根据排泄边界实测的地下水等水位线图确定计算断面位置，并确定各计算断面长度、水力坡降、含水层厚度，按达西公式计算：
Q侧排＝K·I·B·H （7-13）
式中：K为含水层的渗透系数，m/d；I为断面处的水力坡度；B为断面宽度，km；H为潜水含水层厚度，m。灌区2005年侧向补给量计算结果见表7-9。
4.人畜饮水、工业城镇地下水开采量（Q用）
对灌区工业用水、生活用水按各县区进行统计，用水标准按照工业及生活用水定额。近10年灌区人畜饮水、工业城镇地下水开采量见表7-10。
表7-9 灌区2005年侧向排泄量 Table7-9 Lateral discharge in Jinghui Canal Irrigation District in 2005
表7-10 灌区人畜饮水、工业城镇地下水开采量 Table7-10 Canal system leakage recharge in Jinghui Canal Irrigation District
五、地下水均衡分析
根据以上各年水资源量计算结果，计算出区内近10年内平均的补给量、排泄量见表7-11。
表7-11 灌区近10年平均地下水均衡计算结果 Table7-11 Average ground water equilibriumcal culation in Jinghui Canal Irrigation District（2001-2009）
计算区位于泾、渭河河谷阶地区，分布有泾河一至三级；渭河一至二级阶地，其中二级阶地分布面积占总面积的90%左右。在地下水的开采深度范围内，以第四系全新统冲积层为主，次为上更新统冲积层及风积层，西北部边缘及北部与黄士台源接界的局部地区，有第四系全新统洪积扇分布，面积微小。在上述松散岩层中，主要赋存着空隙潜水，仅泾河一级阶地区，由于地层为二元结构（上细下粗），局部地区微具承压性质。从目前及近期发展看，仍以开采第四系全新统冲积含水岩组的潜水为主，高阶地区还涉及上更新统冲积含水岩组（渭河二级、泾河三级阶地区），但此区面积甚小。含水层岩性：在泾河一级阶地区以粗砂、砂砾卵石为主，上覆亚砂士、亚粘士互层，属强富水区，单位涌水量一般为720～168m3/d，由西向东呈递减趋势。渭河一级阶地区以细砂、中粗砂为主，砂砾石次之，间夹数层亚砂士、亚粘士，属强富水区，单位涌水量为360～240m3/d。
泾河二级阶地区，西部及中部为粉细砂，底部有砂砾石分布；东部以粉细砂为主，局部含砂砾石，间夹多层亚粘士、亚砂士，属富水区—中等（偏弱）富水区。由于面积大，受含水层岩相变化及补排条件差异的影响，富水性指标亦有变化。其中单位涌水量：IA区为120～336m3/d；IIB区为120～192m3/d；IIC区为36～72m3/d；泾河三级阶地区以砂层、砂砾石层为主，次为亚粘士、亚砂士夹砂互层，上覆马兰黄士，属中等富水区，单位涌水量为84～120m3/d。渭河二级阶地区为砂砾石、砂层为主，上覆马兰黄士，属中等富水区，单位涌水量为120～240m3/d。
灌区近年在灌溉、开采条件下，浅层地下水的补给源是以灌溉回归水及大气降水的垂向入渗补给为主，约占年总补给量的80%以上（其中灌溉水垂向渗入补给量占年总补给量的38%；降水入渗补给量占年总补给量的55.6%）。其次是径流补给和局部沿河地区的河水侧渗补给，与垂向相比，补给量甚微。潜水的排泄途径也是以垂向为主，径向次之。开采排泄量约占年总排泄量的90%左右，其中人工开采量占89.6%，水平方向排泄量仅占10%左右。灌区绝大部分地区潜水水质较好，唯在灌区北部边缘及高陵以东部分地区，水质较差。灌溉实践中必须注意采用适宜的灌水方式，如渠井汇流、渠井间灌等，以防止有害盐分的集聚。
采用水均衡法对灌区地下水进行评价，各均衡要素计算中所采用的水文地质参数如渗透系数、给水度、导水系数、降水入渗补给系数、干、支渠系渗漏损失系数、灌溉水田间入渗补给系数、井灌回归系数等，是通过非稳定流抽水试验资料以及利用地下水长观资料进行相关分析等方法求得。经过对1951～2009年（灌溉年度）年降水量进行频率计算，选择在灌溉面积、渠灌引水量、地下水开采量、工程设施现状等方面具有代表性及其他资料比较完整的10年（2001年11月1日至2009年10月底），进行全区潜水水均衡计算，计算结果与区内潜水动态变化规律基本一致。说明各种参数的选取及实际资料的采用比较可靠和符合实际。水均衡计算结果：近年多年平均总补给量为2.6767×108m3/a，平均开采地下水资源量为1.6139×108m3/a。
渠系渗漏补给量系指灌溉渠道水位高于地下水位时，渠道水补给地下水的水量。渠系渗漏补给量一般只计算到干、支、斗三级渠道。常见的计算方法有地下水稳定流计算法、经验公式法、渠系入渗补给系数法等。
1.地下水稳定流计算法。计算公式同山丘区全区的多年平均年基流量计算。
2.渠系入渗补给系数法。计算公式为
式中 ?——渠系渗漏补给量，；
?——渠系入渗补给系数；
?——渠首引水量，；
?——渠系渗漏补给地下水系数；
?——渠系有效利用系数。
