一、基础工作
（一）模型识别时段
根据收集到的现有的研究区地下水观测资料，我们取模型识别的时段为2003年7月～2005年6月，取时间步长Δt=1个月，共计24个时间步长。
（二）单元剖分与结点设置
本模型的结点设置和单元剖分遵循下列原则：
1）用于拟合的观测孔需设置结点，角点上设有观测孔的辅助三角形内不得有开采井；用于拟合泉流量的泉眼（点泉）处需设置结点，且该点周围的诸辅助三角单元内不得有抽水井和拟合用的观测孔；泉沟（线泉）应连续地设置结点；
2）各单元尽量接近等边三角形，任一顶角不得大于90°，大小单元间要渐变；
3）河流、渠系的边线尽可能与三角形边线一致，内部的水库水边线及非均质分界线尽可能与差分网格边线一致；
4）单元、结点总数要适当，应与基础资料占有程度及精度相匹配；重点地区应细剖分，缺资料的外围地区可粗剖分。
按照上述原则及研究区的基础资料，将模拟区范围在平面上（即每一个模拟层）剖分为1804个结点，3433个三角单元（图5-12）。模型垂向上共分6个模拟层，6层模拟层的总结点数为1804×6=10824个，总单元数为3433×6=20598个。
图5-12 研究区平面剖分图
（三）模型识别的基础资料
根据研究区的水文地质条件及本章所建立的数学模型确定待定水文地质参数。为此，需要输入计算机的基础数据如下。
1.结点处各模拟层的顶底面标高值
首先对于地面高程，根据已有的研究区的地面标高已知点，并结合MapGIS等高线图，利用Kring插值法生成模型各个结点的地面高程值，使所得到的各结点地面标高值的等值线图（图5-13）与原MapGIS等高线图大致吻合。其次依据本区钻孔（图5-14中圆点）及其分层资料可以得到钻孔在各层的顶底面标高，进而求得各层的厚度，再通过Kring插值得出各结点的厚度。最后根据地面标高及厚度计算出各结点各模拟层的顶底面的标高值。
图5-13 研究区地面标高等值线图
图5-14 研究区钻孔平面分布图
在上述过程中，遇到的一些问题的处理：
1）研究区内大多数钻孔都没有打到下伏基岩，因此其深部的分层需依据“第四系等厚线图”并以部分较完整钻孔的分层资料为标准外推；
2）如图5-14，研究区内钻孔分布不均匀，北部北山区无钻孔分布，对此我们根据研究区最靠近北山的一排钻孔的各层底面标高值，把模型各层底面水平向北推移，从而求得北山各结点的各层厚度；
3）研究区南部的山前地带（洪积扇区）为单一潜水含水层，各层岩性均为砂砾卵石，对此我们参照北部多层含水层的分层标准将其也划分成6层。
由上述方法获得的各模拟层厚度等值线图，见图5-15。
2.初始水头值
本模型在垂向上有6层模拟层，依地下水数值模拟的要求，必须给出各结点每个模拟层的初始水头值。本研究采用“参数-初始水头迭代法”来确定各层的初始水头分布。
3.一类边界水头
区内的一类边界主要是疏勒河南阳镇-双塔水库段，该段疏勒河所在层位各结点的水位值根据实测的河水位值进行插值获得。模型的西边界在模型识别中也作为第一类边界（模型预测时拟改为第二类边界，由模型识别结果给定流量），其水头值根据等水位线外推，并在模型识别中加以校正。另外泉口标高也作为第一类已知水头边界条件输入。
图5-15 各模拟层厚度等值线图
4.第二类边界流量
模型的北边界、东边界、南边界和下边界均作为第二类边界，其中北边界、东边界和下边界为零通量边界；南边界为流量边界，该边界的地下水侧向补给量的计算方法：根据1∶20万水文地质普查获得山区地下水径流模数，再经区域水均衡给定各断面地下水侧向补给量的初始估计值，并在模型识别中校正。
5.降雨入渗系数及降雨量
降水入渗系数首先按照包气带的岩性进行分区，本研究区的入渗系数可分为4个区：①北山区；②北截山区；③粘性土区；④砂砾卵石（砂砾石）区。其中在第三区和第四区降雨入渗系数与水位埋深有关。我们先给定水位埋深为1m，3m，5m和10m处的降雨入渗系数初始值（水位埋深为0m时降雨入渗系数为1.0；水位埋深大于10m时取10m处的降雨入渗系数值），在程序中计算各结点的水位埋深值，并判断其所落区间，再根据线性插值求得该结点的降雨入渗系数，然后经过不断的调试、优选，从而求得相匹配的系数。
本研究降水入渗补给量的计算采用“滞后补给权系数法”。遗憾的是我们未能收集到研究区内三个气象站（玉门镇站、安西站和敦煌站）从2004年1月到2005年6月的月降雨量资料（模拟时段为2003年7月到2005年6月，降雨量资料因考虑滞后补给需提前19个月），因此我们把三个气象站2000～2003年的各月平均降水量数据对应作为2004年1月到2005年6月的各月降雨量资料，再利用所得的三个气象站2001年12月到2005年6月的各月降雨量资料采用平面插值的方法，求得各结点的月降雨量数据。
6.蒸发系数及蒸发量
关于本研究的潜水蒸发规律公式，我们采用的是《甘肃省黑河干流中游平原区包气带水分运移及均衡要素研究》（甘肃二水，1995）一文中得出的经验公式（指数公式）：
河西走廊疏勒河流域地下水资源合理开发利用调查评价
式中：ε为潜水蒸发强度；ε0为水面蒸发强度；h为潜水埋深；b为经验系数（与土质有关），由模型识别确定。
图5-16给出了蒸发皿蒸发与水位埋深的关系，由图中可知，蒸发皿蒸发量在水位埋深0～1m间变化相当大，在此区间地下水埋深减小0.5m，蒸发皿蒸发量可增2～3倍，因此模型识别过程中地面高程与地下水位将直接影响潜水蒸发量的计算。
图5-17为从1997年1月到1999年12月三个气象站（玉门镇站、安西站和敦煌站）的月水面蒸发量柱状图，从图中可知，蒸发量大的月份在每年的4～9月；蒸发量小的月份则在每年的12月至次年2月。需说明的是，本区无土层水的融冻资料。在融冻季节（一般为11月中旬至翌年的3月份），上部种种入渗补给和蒸发被冻结层所阻隔，而冻结、解冻对潜水起上部补排的控制作用。可见此资料的缺乏将影响模拟的精度。
图5-16 潜水蒸发量与水位埋深之间的关系曲线
图5-17 研究区三个气象站1997～1999年各月蒸发量柱状图
7.河流、渠系及田间灌溉水的入渗补给
对于区内河流的入渗补给段，主要有昌马洪积扇顶部三道沟—八道沟河、疏勒河安西段、党河、榆林河和石油河等（图5-18），各入渗补给段的入渗量均按线状量输入。
对于区内渠系入渗补给的处理：各总干、干渠的入渗量按线状量输入（图5-19）；各支、斗渠的入渗量计入各灌区田间灌溉水的入渗量，按面状补给量处理。研究区内有7个灌区（图5-20）：①党河灌区；②西湖灌区；③双塔灌区；④榆林灌区；⑤桥子灌区；⑥昌马灌区；⑦花海灌区。各灌区田间灌溉水的入渗量按单位面积平均分配。
图5-18 研究区水系分布图
8.地下水开采量
区内开采井（包括农灌井和饮水井）很多，且分布比较集中，其开采量难以一一统计，只能通过划分开采区按面井处理。研究区可划分出7个开采区（图5-21）：①党河开采区；②西湖开采区；③双塔开采区；④榆林开采区；⑤桥子开采区；⑥昌马开采区；⑦花海开采区。各开采区开采量的调查是以自然村为单位开展的。
9.点泉和泉沟流量的动态资料
在昌马洪积扇前缘细土平原带，地下水大量以泉的形式排泄。对有系统观测资料（2003年7月～2005年6月）的典型点泉（泉A、泉B和泉C），我们以第一类边界条件（泉口标高即为该点水位标高）模拟其泉流量，然后与实际观测量作拟合。区内泉沟（共22条）溢出量的计算是本模型的一大难点，我们把泉沟线上各结点分别按典型泉的模拟方法来确定其各自泉点的流量，这些泉点形成的总流量即为泉沟的溢出量，也是模型拟合的对象。因我们只收集到2004年7月各泉沟的实测流量资料，故我们只对该月数据进行拟合。各点泉和泉沟的平面位置见图5-22。
10.观测孔地下水头动态
目前，我们收集到的本区系统观测地下水头的动态资料是从2003年7月到2005年6月，共24个月。区内共有53个观测孔用于水头拟合，其中有39个单层观测孔以及14个混合观测孔。观测孔的平面分布如图5-23所示，玉门-踏实盆地内有25个观测孔，安西-敦煌盆地有22个观测孔，花海盆地内有6个观测孔。
图5-19 研究区渠系分布图
图5-20 研究区灌区分布图
图5-21 研究区开采区分布图
图5-22 研究区点泉和泉沟平面分布图
二、水文地质参数的确定
（一）待求的水文地质参数
根据本章所建立的数学模型，并考虑到已有的资料（大气降水量、水面蒸发量以及河流、渠系、田间灌溉入渗补给量、点泉和泉沟的流量和地下水开采量等），确定待求的参数有：各含水层和弱透水层的水平渗透系数、垂直渗透系数、单位储水系数，无压层的重力给水度，降水入渗补给系数，有关潜水蒸发的系数，及田间灌溉入渗补给系数等，而每一参数又可分为若干区。
图5-23 研究区观测孔平面分布图
（二）水文地质参数的确定步骤
在上述数学模型、单元剖分及所选定有关资料的基础上，先按地貌单元及其沉积物类型、岩性特征对待求参数作初步分区，并根据岩性、含水层性质、等水头线、水头动态和抽水试验等资料初步确定各参数的上下限作为约束条件。
我们将前述53个观测孔的水位、点泉和泉沟流量以及自流井流量的动态作为拟合对象，来确定水文地质参数。
在正确建立水文地质概念模型（包括分区）、数学模型及采用正确的仿真数值模拟技术的基础上，观测孔水头拟合的程度反映了所求参数的精度，为此设目标函数：
河西走廊疏勒河流域地下水资源合理开发利用调查评价
式中：p1，p2，…，pn为待求的各水文地质参数；Ng，Nch为对比用的观测孔数和对比时段数；ωh为水头权因子；H（i，j）为模型模拟的i号观测孔j时段的水头，m；Hg（i，j）为i号观测孔j时段的实测水头，m；Nsp，Nat，Nch分别为对比用的泉数、自流井数和对比时段数；ωq为流量权因子；Q（i，j）为模型模拟的i号泉或自流井j时段的流量，m3/d；Qsp/at（i，j）为i号泉或自流井j时段的实测流量，m3/d；λ为观测孔水头拟合所占的比重；（1-λ）为泉和自流井的流量拟合所占的比重。
在各参数上下限约束条件下，求取目标函数G最小，若此时的水位拟合差未达到足够小，则进一步分析其原因，调整相关的分区及参数值，再次求取目标函数最小。如此反复计算、分析，直至求得满意的结果为止。
整个求参步骤见图5-24。
（三）参数的确定
首先粗略地按地貌、沉积类型和岩性等划分各参数分区，给出各参数初始估计值，经过多次调试、优选，最终使各观测孔（不包括2002观测孔）、泉流量及自流井流量的模拟值与实测值达到最佳拟合。优选所得的本区主要水文地质参数见表5-1，与它对应的各参数分区情况见图5-25和图5-26。
图5-24 确定水文地质参数的计算流程图
调试、优选后的降雨入渗补给系数α在北山、北截山区分别为0.03和0.02。粘性土区（第3区）和砂砾卵石（砂砾石）区（第4区）的降雨入渗系数α与水位埋深（h）有关，见图5-26A（其中3代表粘性土，4表砂砾卵石或砂砾石；图中当水位埋深大于或等于10m时，粘性土的降雨入渗补给系数为0.003，砂砾卵石或砂砾石的降雨入渗补给系数为0.001）。
表5-1 分区参数值列表
图5-25 研究区参数分区图
图5-26A 降雨入渗系数与水位埋深间的关系曲线
图5-26B 降雨入渗系数分区图
（四）观测孔的拟合曲线和动态分析
最终得到的各观测孔处模拟水位和实测水位的拟合曲线如图5-27所示。其中，实心圆点表示各观测孔水位的实测值；空心圆点表示模拟水头值；横坐标代表年月；各观测孔孔号见曲线上方（NZKF1表示农综开发1号）。
河西走廊疏勒河流域地下水资源合理开发利用调查评价
河西走廊疏勒河流域地下水资源合理开发利用调查评价
河西走廊疏勒河流域地下水资源合理开发利用调查评价
河西走廊疏勒河流域地下水资源合理开发利用调查评价
河西走廊疏勒河流域地下水资源合理开发利用调查评价
河西走廊疏勒河流域地下水资源合理开发利用调查评价
图5-27 观测孔水头拟合曲线图
为了说明研究区内观测孔水头的拟合情况，将拟合水头的误差统计于表5-2中，用作统计的数据是各拟合观测孔模拟水头与实测水头之差的绝对值（不包括2002观测孔）。
表5-2 观测孔模拟水头与实测水头绝对误差统计表
误差统计表明，ΔH≤0.5m的占总对比数的70.90%，ΔH≤1.0m的占93.57%，总的来说，观测孔水头的拟合效果是较好的，某些观测孔受随机因素（例如附近抽水）的影响拟合误差较大。限于区内地下水开采量和田间灌溉水的入渗量难以统计到各月的数据，而采用各月平均分配的办法，因此大多数受灌溉和开采因素主要影响的观测孔模拟水位年变化比较均匀，波幅较小，没有实测值波动大。
需要说明的是，位于昌马戈壁中部的S97-2及2002号观测孔水位的动态变化，由于昌马水库从2002年底开始蓄水，昌马洪积扇三道沟至八道沟河的入渗量相应减少，导致在模拟时段内，两孔水位一直处于下降趋势（图5-28）。对这两个观测孔水位的拟合，是模型识别中的一个难点。
图5-28 S97-2和2002观测孔水位动态实测曲线
本区以前所做的两个模型（水环所，2000）和（清华大学，2004），都未给出这两个观测孔水位的拟合情况，而S97-2及2002号观测孔水位的拟合，对昌马洪积扇顶部的水文地质分区的划分以及其水文地质参数的确定起着决定性作用。本研究设计了多套方案对此两孔水位进行拟合，但结果表明，若要使2002号观测孔模拟水位曲线与实测水位曲线基本拟合，则2002孔所在水文地质分区的渗透系数值要在500m/d左右。由于该孔岩性未作级配分析及正规抽水试验，我们难以确认该参数值的可取性，根据甘肃省地调院提供的昌马洪积扇砂砾卵石渗透系数的最大值（160m/d左右），我们暂不采用该参数值。究其原因，可能是由于2002号观测孔于20世纪60年代开始使用，至今40多年来一直未洗孔，其观测数据的准确性尚有待确认。故本次研究暂不拟合2002观测孔的水位，只拟合S97-2观测孔，同时确保2002观测孔所在结点的水位值在模型识别时段内呈下降趋势，同时将实测值与模拟值表示在图5-27上。
（五）点泉、泉沟和自流井的流量拟合曲线和动态分析
研究区内各泉点（泉A、泉B和泉C）的流量拟合曲线，如下图5-29所示。其中，实心圆点表示各泉点的实测流量值；空心圆点表示模拟流量值；横坐标代表年月；各泉点标号见曲线上方。
图5-29 泉（泉A、泉B和泉C）流量拟合曲线
由上述拟合曲线可看出，泉B实测流量的数据似乎失真，模拟值与实测值不宜比较；从泉A和泉C的流量拟合曲线可看出，总体上模拟的泉流量变化趋势基本与实测流量相符，但模拟值的波动性没有实测的大，这可能是由于抽水量等源汇项月输入数据精度不足所致。
由于我们只有2004年7月各泉沟的实测流量资料，故只能对该月各泉沟的流量数据进行拟合，拟合曲线如图5-30所示，该月各泉沟流量与实测流量拟合较好。
河西走廊疏勒河流域地下水资源合理开发利用调查评价
图5-30 泉沟流量拟合曲线
图5-31 自流井流量拟合曲线
由上图可看出，模拟所得的各泉沟流量值在每年的6～9月偏小，这是因为研究区内各泉沟均位于玉门-踏实盆地内昌马灌区、桥子灌区和榆林灌区内，而这些灌区的地下水开采均集中在每年的6～9月，从而引起在该时段各灌区的地下水位普遍降低，进而导致各泉沟流量减小。从图5-30还可看出，模拟所得的各泉沟流量基本上是逐年衰减的，这是由于各灌区地下水开采过量引起该区地下水位逐年下降从而导致各泉沟流量逐年衰减。
区内仅一个自流井，其流量的拟合曲线示于图5-31。从图可看出，模拟的自流井流量变化趋势基本与实测流量相符。
（六）模拟时段的潜水面等值线
由于只有2004年6月的实测潜水位等值线图，故取该时段的模拟潜水位等值线（采用克里格插值法绘制）与实测潜水位等值线进行比较（图5-32）。从图中可看出，模拟潜水位等值线与实测潜水位等值线吻合较好，只是在昌马洪积扇顶部，模拟与实测潜水位等值线有所偏差，这是因为洪积扇顶部观测孔较少且该处邻近模型边界，实测值的插值和模拟值都难以精确反映该处地下水流场；北截山、北山及北山山前戈壁均被概化入我们模型内，而实测潜水位等值线图在这些位置均没有绘出潜水位等值线，因此这些地方的模拟潜水位等值线无法与实测潜水位等值线进行比较。在北截山和北山山前地带，模拟与实测潜水位等值线有所偏差，这是由于这些地段缺乏拟合用的观测孔。
由图还可看出，玉门-踏实盆地的地下水由南向北径流，大致在七道沟—饮马农场一带存在地下分水岭，向东经黄花农场流入花海盆地，向西经桥子、踏实向北截山南麓汇集，在昌马洪积扇前缘地下径流以泉的形式溢出地表汇入疏勒河，再经双塔水库流入下游安西-敦煌盆地；安西盆地地下水由东向西径流，在西湖农场—玉门关之间汇党河洪积扇由南向北的地下径流，再经玉门关一带向西径流；花海盆地的地下水由南、东、西三个方向向中北部的干海子汇流（以上地名见图5-33）。地下水垂向径流方面，总体来说在南部山前地下水由上向下流动，至北部细土平原区，转化为向上运动为主。
图5-32 2004年6月实测及模拟潜水位等值线图
图5-33 疏勒河流域模拟区位置图
（七）模拟时段的水均衡分析
由于模拟时段中的2003年和2005年都不完整，故此处仅给出2004年的水均衡资料（表5-3A、5-3B、5-3C和5-3D），作一简单分析。
表5-3A 花海盆地模拟区2004年1月1日～2004年12月31日的水均衡一览表
表5-3B 玉门-踏实盆地模拟区2004年1月1日～2004年12月31日的水均衡一览表
表5-3C 安西-敦煌盆地模拟区2004年1月1日～2004年12月31日的水均衡一览表
表5-3D 疏勒河流域模拟区2004年1月1日～2004年12月31日的水均衡一览表
从表5-3D可看出，研究区地下水的补给主要来源于地表水的入渗，河流、渠系及田间灌溉入渗补给量之和约占总补给量的79.28%，相比之下降凝水入渗对本区地下水的补给非常有限（11.52%）；研究区地下水的排泄主要有蒸发、向疏勒河干流排泄、点泉、泉沟和人工开采等，潜水的蒸发蒸腾消耗了大量的水资源，约占总排泄量的73.55%，区内地下水的开采量占总排泄量的12.77%，向疏勒河干流排泄量和点泉、泉沟及自流井溢出量之和占总排泄量的12.07%。
地下水资源评价中利用均衡法计算的2004年工作区地下水均衡状况见表5-4。
表5-4 工作区2004年1月1日～2004年12月31日均衡计算成果一览表
比较两种方法的计算结果可以看出，除蒸发蒸腾量差别达1.1亿m3/a外，其他各共有项非常近似，这从另一个方面验证了模型的精度。
（八）两种方法水均衡结果的讨论
三维流模型的水均衡计算与区域水均衡计算是独立完成的，两者的结果除蒸发量外，其他诸均衡要素比较接近。蒸发量差别稍大的主要原因有：
（1）计算的方法不同
数值模型是将研究区剖分为1804个小均衡区（因每个小均衡区的潜水埋深和包气带岩性不同），且每个月都做一次水均衡，再将12个月的结果累加成年度的水均衡成果。而区域均衡法则受方法自身的限制，只能按水位埋深的区间值（＜1m，1～3m，3～5m）分区，时间上也只能按年度的始末来计算，未考虑年内可能存在越限的可能。两者的精度水平存在明显的差异。
（2）采用的参数不同
本区缺少潜水蒸发量随潜水埋深和包气带岩性变化的试验资料，数值模型是采用张掖平原堡均衡试验场的经验公式。而区域均衡法则只能按水位埋深的区间值（＜1m，1～3m，3～5m）给出参数，只要水位埋深处在同一区间，蒸发强度就一样，然后考虑有无植被乘以系数。因此参数选取的不同必然导致计算结果的差别。
（3）采用的手段不同
潜水埋深是影响潜水蒸发强度最主要的因素，数值模型是依地面标高减去潜水位而得。其中的地面标高是在1∶25万地形图上选点、插值而得各格点的标高值；潜水位则是数值模拟的结果，其精度取决于模型的仿真性。整个计算由计算机完成。而区域均衡法则只能按有限的观测孔潜水位埋深值的插值而得，由于研究区大、观测孔数少且多集中于绿洲区，大面积的潜水位埋深值靠外推法估计，而且潜水位埋深等值线由手工勾画完成，精度有限。
综上分析，可以保守地估算一下，就蒸发而言，依表6-6的无植被的蒸发强度，＜1m的蒸发强度是1～3m的约4倍；1～3m的蒸发强度也是3～5m的约4倍；5～10m的蒸发强度则是3～5m的约60倍；如果0～5m的某埋深区间的埋深误差0.5m，则有1/4的面积（埋深按平均分布计）误划到另一埋深区间，则此项的蒸发量误差可达一倍。显然，区域均衡法所得埋深值难以保证其误差小于0.5m。
就储存量的改变量而言，它是依潜水位的改变量乘以给水度和面积而得（暂先不说弹性储量）。若潜水位或埋深的误差0.2m，给水度按0.1计，则10000km2面积的释放或存储水量为2亿m3，而研究区的面积达20000km2。显然，区域均衡法所得埋深值难以保证其误差小于0.2m。
基于上述分析，本次大区域的模型计算结果是基本可信的。