降水入渗补给系数

2.4.4.1 降水入渗补给系数变化规律认识

地下水不开采处于完全自然状态时，从长期均衡的角度看，非饱和带厚度和地下水的垂向补给条件基本不变，因而降水入渗补给系数也基本不变。随着地下水开采，潜水和微承压水水位下降，非饱和带厚度增大，地下水的垂向补给条件发生变化，降水入渗补给系数也随之发生变化。因此，分析和研究随着地下水开采程度的增长，降水入渗补给系数的变化规律，对于充分认识潜水和潜水-微承压水的垂向补给变化，充分开发利用潜水和潜水-微承压水，科学地控制地下水的开采，以达到持续利用和保护环境的目的是十分重要的。

尤其在我国北方地区，地下水开采程度比较高，地下水开采经历了一个长时间发展过程，积累了大量地下水动态等相关资料，为这方面的调查研究创造了有利条件。

分析研究的重要内容:

(1)年降水量P年、次降水量Px与降水入渗补给量和降水入渗补给系数的相互关系，绘制α-P年、α-Px关系曲线。分区选择历年不同降水量资料及其附近地下水观测点的水位升幅值，求取降水入渗补给系数。并绘制相关曲线图。(注意非饱和带岩性的相似性)。

(2)地下水埋深(h)与降水入渗补给系数的相互关系。

(3)不同非饱和带岩性条件下，降水入渗补给系数与地下水埋深的相互关系，绘制其相关曲线。如图2.4.1 所示。

图2.4.1 不同岩性降水入渗补给系数α-地下水埋深(h)相关曲线

2.4.4.2 降水入渗补给系数评估计算方法

在长期地下水调查研究实践中已积累了许多方法，但各种方法也有其局限性，实际工作中要了解这

些方法的特点、存在问题和相互差异，结合各地区自然条件和已有资料来选择合适的计算方法。

2.4.4.2.1 利用地下水长观动态资料求降水入渗补给系数由于地下水长期观测区在区域上数量多、分布较均匀，各类代表性地区大都没有观测点，因此，用地

下水长观动态资料求取降水入渗补给系数是普遍且效果较好的方法。其计算方法主要有两种:

(1)年水位升幅累积法，计算年降水入渗补给系数

这个方法的前提是一些平原区地下水侧向流动较缓慢，天然条件下，地下水位升幅完全代表了地下水含水层所获得的降水入渗补给。因此年降水入渗补给系数为降水所引起的地下水升幅之和乘以给水度被年降水量除。

地下水资源调查评价技术方法汇编

式中:μ———给水度;

Δhi———降水引起的次水位升幅;

N———全年降水次数，i＜N;

∑Pi=P年———年降水总量;

Ni———年内降水引起水位升幅的有效补给的次数，Ni＜N。

(2)前期影响降水量法

这个方法需要研究程度比较高，尤其适用于非饱和带土壤水运移规律和参数研究比较有基础的地区。

方法要求仔细研究历次降水过程和补给，以及前次降水的影响。首先计算次降水入渗量，次降水入渗补给系数，再换算成年降水入渗补给系数。计算公式:

地下水资源调查评价技术方法汇编

式中:Pai———前期降水影响量(它概括反映了无效降水量和非饱和带土层含水量对次降水入渗补给量的影响);

k———影响系数，取值0.85～0.95，一般取平均值0.9;

Pi———本次降水量;

Pi-1———前次降水量;

t———距本次降水的天数，可由本次降水向前推15～20天。次降水入渗补给系数:

地下水资源调查评价技术方法汇编

式中:μ———给水度;

Δhi———由本次降水入渗补给形成的水位升幅;

Pi———本次降水量;

Pai———本次降水的前期降水影响量。

需要计算年降水入渗补给系数时，须将次降水入渗补给系数换算成年降水入渗补给系数。计算公式如下:

地下水资源调查评价技术方法汇编

利用动态资料求取降水入渗补给系数注意的问题:

·由于地下水在含水层中多年循环和调节补给的结果，地下水资源不完全对应于每年的补给量，而是在一个有代表性的气象周期内平均值的概念。需要进行多年调节计算时，要采用相应频率的降水量数据。年降水量小于400mm以下，降水入渗补给系数将明显减少。

·计算时注意分离出非降水因素引起的水位升幅值，如由于河水灌溉引起的水位上升影响等。

·地下水长期观测一般为5日一次，也有10日一次的情况。在雨季，这样的间隔对降水入渗补给系数计算影响较大。如果有试验区的日观测资料，可以用日观测资料计算后对5日、10日观测资料计算结果加以修正。

2.4.4.2.2 均衡试验场观测试验求取降水入渗补给系数

国土资源部和其他相关行业部门都曾在各地建设了一些试验场，这些试验场目的不尽相同，但不少试验场都设置有气象和降水入渗观测。收集这些资料对调查地下水补给和确定参数极为重要。

地下水均衡研究中，对降水入渗补给和蒸发的观测主要有三种方法:第一种是马里奥特瓶法，第二种是称重法，这两种方法都将降水入渗与蒸发分别观测计算。因此，马氏瓶中接收的水量和称重得到的水量都不能反映自然界地下水在含水层中补给-蒸发作用同时进行所获得的真实补给量和形成的水位升幅降幅。因此只有将相同时间段内降水入渗量减去蒸发量，才能得到真正的降水入渗补给量。第三种方法是非饱和带水分运移观测试验，测量降水通过非饱和带下渗补给潜水的水量，通常通过零通量面法和瞬时流量法试验观测计算获得。

以上试验取得的都是点上的年降水入渗系数，用到较大面积的计算分区时，会有差异。

2.4.4.2.3 模型模拟反求降水入渗补给系数

求得的降水入渗补给系数经过模型模拟均衡计算，通过几年均衡量的校核，可信度比较高。要注意的是，模型计算所得参数是模型计算中参数分区面上的平均降水入渗补给系数值。开采量数据的准确性对反求参数的大小影响极大，应首先考核开采量的准确性和开采量的年内分配合理性。

2.4.4.2.4 岩溶小泉域排泄量反求降水入渗补给系数

对边界比较清楚的全排型岩溶小泉域可以统计泉域内的泉和集中排泄带的泉的流量，除以小泉域内面降水量，求取降水入渗补给系数。

此参数可引入区内其他岩溶分布区应用，须注意的是:

(1)泉域内地下水已开采利用时，对野外测定的泉流量要进行地下水开采量还原。

(2)泉域内第四系覆盖薄，降水易渗入到岩溶含水层中，且有农业灌溉的地区，则需注意计算中应扣除部分灌溉入渗补给量。

2.4.4.2.5 利用观测孔组资料用有限差分法求降水入渗补给系数

当区内长期观测孔数量不足或分布不匀时，可以用若干长观孔，或者调查时临时选定若干计算时段，计算孔组，组织短期观测试验后用有限差分法计算降水入渗补给系数。

为减少干扰，提高计算精度，计算时段的选择极为重要。一般选择无开采、无灌溉的时段进行为宜。

以一个观测孔为中心，与周边若干观测孔可以连成几个三角形，以各三角形的每一边中点所作垂线，交点连成的多边形是观测孔的均衡区，选定Δt计算时段内单位面积垂向补给量W(图2.4.2)。

图2.4.2 观测孔均衡区划分

地下水资源调查评价技术方法汇编

式中:μ———给水度;

Δhi———中心孔在Δt计算段内水位上升值;

T———导水系数;

Δt———计算时段;

A1———中心孔均衡面积;

———任一观测孔平均水位高程;

h1———中心孔平均水位高程;

L1i———中心孔与周边任一观测孔i的距离;

l1i———中心孔与周边各观测孔连线中点垂线组成的泰森多边形的边长。

用此方法也可计算区域平均蒸发量。

当计算期内有开采井时，上述计算公式右端应加一项:

地下水资源调查评价技术方法汇编

式中:m———开采井个数;

Qp———任一井抽水量;

A1———中心孔均衡面积;

aβp———P井与中心孔对边的距离;

Lβ———中心孔与对边的距离;

———流量分配系数，

水资源数值模型校验及现状地下水均衡分析

c feng.for五层准三维流（潜/弱/承/弱/承，对弱透水层只考虑越流不考虑弹性释放），调参用

c 使潜、承、承上下节点、上下单元完全重合，且总单元数ms＝3×ms1，总节点数ids＝3×ids1

c 潜单元编号由1-ms1，中承单元编号由ms1+1-ms2，深承单元编号由ms2+1-ms

c 10-13行：可调数组的各个参数预先赋值，不同计算区只要改变这些参数即可用此程序，其中潜单元个数＝ms1个，中承单元个数＝ms2-ms1个，深承单元个数＝ms-ms2个；潜节点编号由1-ids1，中承节点编号由ids1+1ids2，深承节点编号由ids2+1-ids，n：未知水头节点数；mi：开采节点数；igs：拟合点总数；igs1：潜拟合点数，igs2igs1：中承拟合点数，icqs：参数区个数：ihv：计算时段数，jbn：纯外引地表水灌溉节点个数；nn1：潜已知水头节点个数；nn：潜+中+深已知水头节点总数

parameter（ms1＝108，ids1＝67，n1＝59，mi1＝17，igs1＝11，icqs1＝8）

parameter（ms2＝216，ids2＝134，n2＝126，mi2＝47，igs2＝25，icqs2＝16）

parameter（ms＝324，ids＝201，n＝193，mi＝69，igs＝35，icqs＝24）

parameter（ihv＝24，jbn＝9，nn1＝8，nn＝24）

c h0：时段初水头（m）；hed：时段末水头（m）；in：各单元三节点编号（必须按小中大顺序排）

c fh：拟合点各时段计算水位（m）；sh：拟合点各时段实测水位（m）

dimension h0（ids），hed（ids），in（ms，3），fh（igs，ihv），sh（igs，ihv）

c idyh：导水矩阵工作单元；zb：各节点x，y坐标（输入坐标为图面mm数）；igdh：拟合点节点号；q：开采井各时段水量（抽为正，注为负，流出边界为正，流入为负）（m3/d）

dimension idyh（ids，9），zb（ids，2），igdh（igs），q（mi，ihv）

c d：各节点导水矩阵；ifdh：开采井节点号；s：拟合点拟合误差（m）；mqh：各单元的参数区号；e：各节点释（储）水矩阵dimension ifdh（mi），d（ids，9），s（igs，ihv），mqh（ms），e（ids）

c fqc：各参数区最多4个参数之值；sss：各节点水头总变幅，由t0时刻起算，（m）；h：初始流场，gh潜水层独有的灌溉入渗系数

dimension fqc（icqs，4），sss（ids），h（ids），gh（icqs1）

c yo：各节点越流矩阵工作单元；jbiao：纯引地表水灌溉的节点编号，共有jbn个节点；qxia：灌溉抽取地下水均铺在灌溉区的水层厚度；gq：灌溉抽出的地下水水层厚度；gq＝qxia；w：各时段潜水各节点的降水灌溉回归水量

dimension yo（ids），jbiao（jbn），qxia（ihv），gq（ihv），w（ids1，ihv）

c zzz：潜水各节点底板标高（m）；ep0：潜水各节点含水层厚（m）；ep1：潜水各单元含水层厚均值（m）；x：各时段降水量（m/时段）

dimension zzz（ids1），ep0（ids1），ep1（ms1），x（ihv）

c tl：各时段步长值（天）；bc：三角单元几何量（bi，ci，bj，cj，bk，ck）

dimension tl（ihv），bc（3，2）

c 分别是潜（h0nn1），中（h0nn2），深（h0nn3）层已知变水头节点各时段的水头

dimension h0nn1（nn1，ihv），h0nn2（nn1，ihv），h0nn3（nn1，ihv）

c file1（igs）*8：定义8字符的字符串igs个；filee*8：定义8字符的文件名

character file1（igs）*8，file2（igs）*8，filee*8

write（*，23）

23 format（1x，’zz1＝？zz2＝？迭代因子：zz1潜水0—1，zz2承压水1—2’）

c zz1：潜水为亚松弛系数，一般取0.85为好；zz2：承压水为超松弛系数，当ids＜300时取1.2-1.3，当ids＞300时取1.3-1.5为好

zz1＝0.85

zz2＝1.25

write（*，24）ms1，ids1，n1，mi1，igs1，icqs1

write（*，24）ms2，ids2，n2，mi2，igs2，icqs2，ihv

write（*，24）ms，ids，n，mi，igs，icqs

24 format（5x，12i5）

open（1，file＝’fqc’，status＝’old’）

c 首先从’fqc’中读取潜水的区号m10（空读），参数1（K），参数2（μ），参数3潜-中的越流系数（k’/m’），参数4降水入渗系数（α），gh灌溉入渗系数

read（1，*）（m10，fqc（i，1），fqc（i，2），fqc（i，3），fqc（i，4），gh（i），

* i＝1，icqs1）

c 接着从’fqc’中读取中层水的区号m10（空读），参数1（T），参数2（μ*）

read（1，*）（m10，fqc（i，1），fqc（i，2），i＝icqs1+1，icqs2）

c 接着从’fqc’中读取深层水的区号m10（空读），参数1（T），参数2（μ*），参数3：中-深的越流系数（k’/m’）

read（1，*）（m10，fqc（i，1），fqc（i，2），fqc（i，3），i＝icqs2+1，icqs）

close（1）

c 参数3：中-深的越流系数（k’/m’）同样用于中层水

do 54 i＝1，icqs1

fqc（icqs1+i，3）＝fqc（i，3）

54 continue

c qxia：灌溉抽取地下水平均铺在灌溉区的水层厚度

open（1，file＝’gq’，status＝’old’）

read（1，*）（qxia（i），i＝1，ihv）

close（1）

c gq：灌溉抽出的地下水水层厚度

open（1，file＝’gq’，status＝’old’）

read（1，*）（gq（i），i＝1，ihv）

close（1）

c jbiao：纯引地表水灌溉的节点编号，共有jbn个节点

open（1，file＝’jbiao’，status＝’old’）

read（1，*）（jbiao（i），i＝1，jbn）

close（1）

c’in’：各单元三节点编号文件（必须按小中大顺序排），ia为单元号空读

open（1，file＝’in’，status＝’old’）

read（1，*）（ia，（in（i，j），j＝1，3），i＝1，ms）

close（1）

c igdh：拟合点节点编号

open（1，file＝’igdh’，status＝’old’）

read（1，*）（igdh（i），i＝1，igs）

close（1）

c mqh：各单元的参数区号，ia为单元号空读

open（1，file＝’mqh’，status＝’old’）

read（1，*）（ia，mqh（i），i＝1，ms）

close（1）

c ifdh：开采井节点号

open（1，file＝’ifdh’，status＝’old’）

read（1，*）（ifdh（i），i＝1，mi）

close（1）

c tl（i）：各时段值（天/时段，每月为1个时段）

open（1，file＝’tl1’，status＝’old’）

read（1，*）（tl（i），i＝1，ihv）

close（1）

c 每个时段的大气降水量（m），x（i）/dt则换算为（m/d）

open（1，file＝’x’，status＝’old’）

read（1，*）（x（i），i＝1，ihv）

close（1）

c’zb’：各节点x，y坐标文件（此处为1：2.5万图面的mm数），ia为节点号空读

open（1，file＝’zb’，status＝’old’）

read（1，*）（ia，（zb（i，j），j＝1，2），i＝1，ids）

close（1）

c 把各节点x，y坐标1：2.5万图面的mm数换算为实地的m数

do 888 i＝1，ids

do 888 j＝1，2

zb（i，j）＝zb（i，j）*25

888 continue

c 潜水开采井各时段的开采量（m3/d），空读ia开采井个数，空读ib开采井所在节点编号

open（1，file＝’qq’，status＝’old’）

read（1，*）ia，（（ib，q（i，j），i＝1，mi1），j＝1，ihv）

close（1）

c 中层水开采井第一时段的开采量（m3/d），空读ia开采井个数，空读ib开采井所在节点编号

open（1，file＝’qz1’，status＝’old’）

read（1，*）ia，（ib，q（i，1），i＝mi1+1，mi2）

close（1）

c 中层水开采井第13时段的开采量（m3/d），空读ia开采井个数，空读ib开采井所在节点编号

open（1，file＝’qz2’，status＝’old’）

read（1，*）ia，（ib，q（i，13），i＝mi1+1，mi2）

close（1）

c 深层水开采井第一时段的开采量（m3/d），空读ia开采井个数，空读ib开采井所在节点编号

open（1，file＝’qs1’，status＝’old’）

read（1，*）ia，（ib，q（i，1），i＝mi2+1，mi）

close（1）

c 深层水开采井第13时段的开采量（m3/d），空读ia开采井个数，空读ib开采井所在节点编号

open（1，file＝’qs2’，status＝’old’）

read（1，*）ia，（ib，q（i，13），i＝mi2+1，mi）

close（1）

c 中、深层水开采井第2-12时段的开采量等于第一时段的开采量

do 887 j＝2，12

do 886 i＝mi1+1，mi

q（i，j）＝q（i，1）

886 continue

887 continue

c 中、深层水开采井第14-ihv时段的开采量等于第13时段的开采量

do 884 j＝14，ihv

do 883 i＝mi1+1，mi

q（i，j）＝q（i，13）

883 continue

884 continue

c sh：拟合点各时段实测水位（m），空读ia拟合点所在节点编号

open（3，file＝’sh’，status＝’old’）

read（3，*）（ia，（sh（i，j），j＝1，ihv），i＝1，igs）

close（3）

c qc：潜水初始流场，空读is节点编号

open（1，file＝’qc’，status＝’old’）

read（1，*）（is，h（i），i＝1，ids1）

close（1）

c zc：中层水初始流场，空读is节点编号

open（1，file＝’zc’，status＝’old’）

read（1，*）（is，h（i），i＝ids1+1，ids2）

close（1）

c sc：深层水初始流场，空读is节点编号

open（1，file＝’sc’，status＝’old’）

read（1，*）（is，h（i），i＝ids2+1，ids）

close（1）

c h0nn：分别是潜（h0nn1），中（h0nn2），深（h0nn3）层已知变水头节点各时段的水头，空读ii节点编号

open（1，file＝’h0nn’）

read（1，*）（ii，（h0nn1（i，ikv），ikv＝1，ihv），i＝1，nn1）

read（1，*）（ii，（h0nn2（i，ikv），ikv＝1，ihv），i＝1，nn1）

read（1，*）（ii，（h0nn3（i，ikv），ikv＝1，ihv），i＝1，nn1）

close（1）

c zzz：潜水各节点底板标高（m），空读ii节点编号

open（1，file＝’zzz’）

read（1，*）（ii，zzz（i），i＝1，ids1）

close（1）

c’file1’存放igs个字符串（拟合节点名称占4字符，及观测井的原编号再占4字符）

open（1，file＝’file1’，status＝’old’）

read（1，110）（file1（i），i＝1，igs）

close（1）

c’file2’存放igs个字符串（拟合节点名称占4字符，加.dat后缀再占4字符）

open（1，file＝’file2’，status＝’old’）

read（1，110）（file2（i），i＝1，igs）

close（1）

110 format（10a8）

c 计算开始前先把全部节点的时段末刻水头hed（i），时段初刻水头h0（i）赋初值h（i），以使开始迭代时hed（i），h0（i）不为零

do 1993 i＝1，ids

hed（i）＝h（i）

1993 h0（i）＝h（i）

c 对中、深层承压水导水矩阵工作单元idyh，释水矩阵e，越流矩阵yo，导水矩阵d，赋0

do 25 i＝1+ids1，ids

do 25 j＝1，9

idyh（i，j）＝0

e（i）＝0.0

yo（i）＝0.0

25 d（i，j）＝0.0

c sum2用来累计计算区总面积（用承压水总面积代替），先赋零

sum2＝0

c 对承压水逐个单元计算几何量及导水、释水、越流矩阵

c 对承压水第ip单元的三个节点号依次赋给i，j，k及i1，j1，k1

do 80 ip＝ms1+1，ms

i＝in（ip，1）

j＝in（ip，2）

k＝in（ip，3）

i1＝i

j1＝j

k1＝k

c idyh（i1，9）存放与i1点同单元的所有节点号，最多9个，可以用不完，即i1点的idyh可以有几个idyh（i1，i2）＝0；当计算第ip单元时，i1点的idyh由i1占第1个（i2＝1）位置，j1，k1只能占i2＝2，3，…，9 中的位置；且先占者排前，193行：使j1，k1分两轮到idyh中找位置；当j1找时，195行发现i1的idyh中已有j1，则跳到j1对i1的96循环体头，换为k1对i1循环；当196行没发现已有j1，且196行判断此位不空时返回96循环体头找下一个位置，当碰到第1个空位时，由j1占据，然后跳到j1对i1循环体头，换为k1对i1循环

do 90 iv＝1，3

idyh（i1，1）＝i1

do 94 iu＝j1，k1，k1-j1

do 96 i2＝2，9

if（idyh（i1，i2）.eq.iu）goto 94

if（idyh（i1，i2）.ne.0）goto 96

if（idyh（i1，i2）.eq.0）idyh（i1，i2）＝iu

goto 94

96 continue

94 continue

c 把ip单元的下一个节点j提为i1，k提为j1，i降为k1，然后返回90循环头处理ip单元的下一个节点，循环3次则ip单元中i，j，k 3个点的idyh全部找完

iuu＝i1

i1＝j1

j1＝k1

k1＝iuu

90 continue

c 第ip单元中i，j，k三节点的x坐标赋给xi，xj，xk，y坐标赋给yi，yj，yk

xi＝zb（i，1）

xj＝zb（j，1）

xk＝zb（k，1）

yi＝zb（i，2）

yj＝zb（j，2）

yk＝zb（k，2）

c 第ip单元三角形面积ss

ss＝abs（（xj-xi）*（yk-yi）-（xk-xi）*（yj-yi））*0.5

c 累加各单元面积（这里的sum2累加之后是中、深层两层面积之和）

sum2＝sum2+ss

c 第ip单元的bi～bc（1，1），ci～bc（1，2），bj～b（2，1），cj～b（2，2），bk～b（3，1），ck～b（3，2）

bc（1，1）＝yj-yk

bc（1，2）＝xj-xk

bc（2，1）＝yk-yi

bc（2，2）＝xk-xi

bc（3，1）＝yi-yj

bc（3，2）＝xi-xj

c 第ip单元所在参数区号赋给jv，参数1（T）赋给txy，参数2（μ）赋给ts，参数3（k’/m’）赋给rmk

jv＝mqh（ip）

txy＝fqc（jv，1）

ts＝fqc（jv，2）

rmk＝fqc（jv，3）

c 第ip单元三节点i，j，k的释水矩阵元素c（ii，ii）及c（ii，jj）……式中没包括1/（2Δt）时，随着ip＝ms1+1，ms的单元循环而对有关单元求其和

e（i）＝-ts*ss/3.0+e（i）

e（j）＝-ts*ss/3.0+e（j）

e（k）＝-ts*ss/3.0+e（k）

c越流矩阵元素，第ip单元1m水头差时的越流量（m3/d/m）均分给三节点i，j，k，随着ip＝ms1+1，ms的单元循环而对有关单元求其和

yo（i）＝yo（i）+rmk*ss/3.0

yo（j）＝yo（j）+rmk*ss/3.0

yo（k）＝yo（k）+rmk*ss/3.0

c 计算ip单元各节点的导水矩阵元素d（i，j），（i＝1，2，3）（j＝1，2，3）

do 100 iu＝1，3

do 104 iuu＝1，3

c 当iu＝1时，即ip单元的i节点，分别计算iuu＝1，2，3，即ip单元的i，j，k节点对i节点的ai，ai即公式**的没求和部分

i＝in（ip，iu）

j＝in（ip，iuu）

ai＝txy*（bc（iu，1）*bc（iuu，1）+bc（iu，2）*bc（iuu，2））/ss/4.0

c 在ip单元的三个节点中，排除第3点，只让第1，2点（i，j点）的ai加入i点的导水矩阵元素d（i，j）中，第243行的j可分别轮到i，j，k三点，但第245句的1～9个中，仅有i，j点在i点的idyh中，此句排除了第3点加入i点的导水矩阵元素d（i，j）中随着ip＝ms1+1，ms的循环，到251句时承压水导水矩阵已完全形成

do 106 k＝1，9

if（j.eq.idyh（i，k））d（i，k）＝ai+d（i，k）

106 continue

104 continue

100 continue

80 continue

c 至此，中、深层几何量计算完，以下开始时段循环；时段循环中把潜水也加进来

c 在时段循环中，潜水的导水矩阵，释水矩阵，越流矩阵与含水层厚度有关，所以在时段循环中完成，另外，承压水与时段有关的部分也在时段循环中完成

do 280 ikv＝1，ihv

write（*，*）’ikv＝’，ikv

c 该时段步长赋给dt

dt＝tl（ikv）

c 以下开始计算潜水末刻流场

c 潜水已知水头节点的已知水头赋给本时段初、末刻

do 32 i＝1，nn1

hed（n1+i）＝h0nn1（i，ikv）

h0（n1+i）＝h0nn1（i，ikv）

32 continue

c 先计算潜水几何量，因为潜水的含水层厚度M 随时段而变，所以放在时段循环内

do 33 i＝1，ids1

h0（i）＝hed（i）

c 本时段潜水各节点的含水层厚度ep0（i）用时段初刻水位＝hed（i）减去底板标高zzz（i）求得

ep0（i）＝hed（i）-zzz（i）

c本时段潜水各节点的灌溉回归水量加降水回归水量w（i，ikv）先赋零，以备求和之用

w（i，ikv）＝0.0

33 continue

c 该时段潜水各单元的含水层厚度用三节点之均值

do 299 i＝1，ms1

299 ep1（i）＝（ep0（in（i，1））+ep0（in（i，2））+ep0（in（i，3）））/3.

do 2050 i＝1，ids1

c 同前，对潜水导水矩阵工作单元idyh，释水矩阵e，越流矩阵yo，导水矩阵d，先赋0

do 2050 j＝1，9

idyh（i，j）＝0

e（i）＝0.0

yo（i）＝0.0

2050 d（i，j）＝0.0

do 8015 ip＝1，ms1

i＝in（ip，1）

j＝in（ip，2）

k＝in（ip，3）

i1＝i

j1＝j

k1＝k

do 9015 iv＝1，3

idyh（i1，1）＝i1

do 9415 iu＝j1，k1，k1-j1

do 9615 i2＝2，9

if（idyh（i1，i2）.eq.iu）goto 9415

if（idyh（i1，i2）.ne.0）goto 9615

if（idyh（i1，i2）.eq.0）idyh（i1，i2）＝iu

goto 9415

9615 continue

9415 continue

iuu＝i1

i1＝j1

j1＝k1

k1＝iuu

9015 continue

xi＝zb（i，1）

xj＝zb（j，1）

xk＝zb（k，1）

yi＝zb（i，2）

yj＝zb（j，2）

yk＝zb（k，2）

ss＝abs（（xj-xi）*（yk-yi）-（xk-xi）*（yj-yi））*0.5

jvv＝mqh（ip）

dyk＝fqc（jvv，4）

dyk1＝gh（jvv）

do 776 ig＝1，3

ngh＝in（ip，ig）

do 777 jb＝1，jbn

if（ngh.eq.jbiao（jb））goto 778

777 continue

w（ngh，ikv）＝w（ngh，ikv）-qxia（ikv）*ss/3

778 w（ngh，ikv）＝w（ngh，ikv）+x（ikv）/dt*dyk*ss/3+gq（ikv）*dyk1*ss/3

776 continue

bc（1，1）＝yj-yk

bc（1，2）＝xj-xk

bc（2，1）＝yk-yi

bc（2，2）＝xk-xi

bc（3，1）＝yi-yj

bc（3，2）＝xi-xj

jv＝mqh（ip）

txy＝fqc（jv，1）*ep1（ip）

ts＝fqc（jv，2）

rmk＝fqc（jv，3）

e（i）＝-ts*ss/3.0+e（i）

e（j）＝-ts*ss/3.0+e（j）

e（k）＝-ts*ss/3.0+e（k）

yo（i）＝yo（i）+rmk*ss/3.0

yo（j）＝yo（j）+rmk*ss/3.0

yo（k）＝yo（k）+rmk*ss/3.0

do 1001 iu＝1，3

do 1041 iuu＝1，3

i＝in（ip，iu）

j＝in（ip，iuu）

ai＝txy*（bc（iu，1）*bc（iuu，1）+bc（iu，2）*bc（iuu，2））/ss/4.0

do 1061 k＝1，9

if（j.eq.idyh（i，k））d（i，k）＝ai+d（i，k）

1061 continue

1041 continue

1001 continue

8015 continue

c 至此，该时段潜水几何量计算完

c 时段末水头迭代计数器ikv2，最大误差记录amax

ikv2＝0

9881 amax＝0.0

ikv2＝ikv2+1

write（*，34）ikv，ikv2

34 format（3x，i4，20x，’ikv2＝’，i4）

c 开采井点录用计数器iq

iq＝1

c 对潜水n1个未知水头节点逐点计算该时段末刻水头

do 4002 i＝1，n1

c i节点常数项res：减该时段i点降水灌溉回归量，减i点该时段来自承压水的越流量（越流量计算采用时段初承压水水头，时段末潜水水位），减i点该时段储存量的增加量

res＝-w（i，ikv）-（hed（i+ids1）-hed（i））*yo（i）-e（i）*（hed（i）-h0（i））/dt

if（i.eq.ifdh（iq））then

c 如果i节点恰是开采节点，则i节点常数项res中再加上该时段i点的开采量，然后开采井点录用计数器iq加1，开采井点只能被1个节点录用

res＝res+q（iq，ikv）

iq＝iq+1

endif

c 把与i节点共单元的k节点的导水矩阵元素（k＝1，2，…，最多可9，当k＝1时即i点自己）依次加进i节点常数项res中，k节点的导水矩阵元素采用上一轮迭代出的hed（k）计算

do 3002 j＝1，9

k＝idyh（i，j）

if（k.eq.0）goto 3002

res＝res+d（i，j）*hed（k）

3002 continue

c 把常数项除以：（i节点对i节点自身导水矩阵元素d（i，1）的负数，加，i节点储水阵元素e（i）/dt）

res＝res/（-d（i，1）+e（i）/dt）

c 把i节点的res乘以亚松弛系数zz1，加给上一轮迭代出的hed（i）中，作为这一轮迭代出的hed（i）：

hed（i）＝res*zz1+hed（i）

c 把1-n1个节点中最大的res挑出，并把其点号记到imax1中：

if（abs（res）.le.amax）goto 4002

amax＝abs（res）

imax1＝i

c 循环返回，继续下一个节点：

4002 continue

988 amax＝0.0

232 continue

write（*，*）’amax1＝’，amax，imax1

if（amax.gt.9.999999e-03）goto 9881

c 至此，潜水本时段末刻流场计算完，以下开始中层水本时段末刻流场的计算：

do 232 i＝1，nn1

hed（n2+i）＝h0nn2（i，ikv）

h0（n2+i）＝h0nn2（i，ikv）

iqq＝mi1+1

do 400 i＝ids1+1，n2

c 其中的越流流入量采用（潜水本时段末水位hed（i-ids1）减承压水本时段末水位hed（i））之差hedd计算

hedd＝hed（i-ids1）-hed（i）

heddd＝hed（i）-hed（i+ids1）

res＝-hedd*yo（i）+heddd*yo（i+ids1）-e（i）*（hed（i）-h0（i））/dt

50 if（i.eq.ifdh（iqq））then

res＝res+q（iqq，ikv）

iqq＝iqq+1

endif

do 300 j＝1，9

k＝idyh（i，j）

if（k.eq.0）goto 300

res＝res+d（i，j）*hed（k）

300 continue

res＝res/（-d（i，1）+e（i）/dt）

hed（i）＝res*zz2+hed（i）

if（abs（res）.le.amax）goto 400

amax＝abs（res）

imax2＝i

400 continue

write（*，*）’amax2＝’，amax，imax2

if（amax.gt.9.999999e-03）goto 988

c 至此，中层水全部节点末刻水头迭代完一轮，节点误差最大者如果 ＞0.01m，则988 句进行下一轮迭代，如果≤0.01m，则迭代结束，中层水本时段末刻流场计算完，以下开始深层水本时段末刻流场的计算：

989 amax＝0

iqqq＝mi2+1

do 332 i＝1，nn1

hed（n+i）＝h0nn3（i，ikv）

h0（n+i）＝h0nn3（i，ikv）

332 continue

do 401 i＝ids2+1，n

hedd＝hed（i-ids1）-hed（i）

res＝-hedd*yo（i）-e（i）*（hed（i）-h0（i））/dt

5011 if（i.eq.ifdh（iqqq））then

res＝res+q（iqqq，ikv）

iqqq＝iqqq+1

endif

do 301 j＝1，9

k＝idyh（i，j）

if（k.eq.0）goto 301

res＝res+d（i，j）*hed（k）

301 continue

res＝res/（-d（i，1）+e（i）/dt）

hed（i）＝res*zz2+hed（i）

if（abs（res）.le.amax）goto 401

amax＝abs（res）

imax3＝i

401 continue

write（*，*）’amax3＝’，amax，imax3

if（amax.gt.9.999999e-03）goto 989

c至此，深层水全部节点本时段末刻水头迭代完一轮，节点误差最大者如果＞0.01m，则返回989句进行下一轮迭代，如果≤0.01m，则迭代结束，深层水本时段末刻流场计算完

c以下开始把本时段末刻水位潜，中，深层水全部节点本时段末刻水头hed赋给下时段初刻水位h0，三层拟合点本时段末刻水位hed存入fh

189 do 281 i＝1，ids

h0（i）＝hed（i）

281 continue

do 190 j＝1，igs

mp＝igdh（j）

fh（j，ikv）＝hed（mp）

190 continue

280 continue

c 至此，三层全包括的时段循环完

c 屏幕输出计算区总面积sum2/2，（m2）；sum2是中、深层两层的累加和：

write（*，*）’sum＝’，sum2/2，’m2’

do 1926 i＝1，ids

1926 sss（i）＝hed（i）-h（i）

c’ok1’输出潜水的最终流场：节点号i、坐标zb、最终流场hed、该节点0-ihv时间（最终流场比初始流场）的水位总变幅：

open（1，file＝’ok1’）

write（1，1927）（i，zb（i，1），zb（i，2），hed（i），sss（i），i＝1，ids1）

close（1）

c’ok2’输出中层水的最终流场：节点号i、坐标zb、最终流场hed、该节点0-ihv时间（最终流场比初始流场）的水位总变幅：

open（1，file＝’ok2’）

write（1，1927）（i，zb（i，1），zb（i，2），hed（i），sss（i），i＝ids1+1，ids2）

close（1）

c’ok3’输出深层水的最终流场：节点号i、坐标zb、最终流场hed、该节点0-ihv时间（最终流场比初始流场）的水位总变幅：

open（1，file＝’ok3’）

write（1，1927）（i，zb（i，1），zb（i，2），hed（i），sss（i），i＝ids2+1，ids）

close（1）

1927 format（i4，2f7.0，2f7.2）

c 计算出拟合点各时段误差（计算水位-实测水位）s（i，iuv）

do 8 iuv＝1，ihv

do 2003 i＝1，igs

s（i，iuv）＝fh（i，iuv）-sh（i，iuv）

2003 continue

8 continue

c 输出拟合点各时段误差

open（1，file＝’gan.dat’）

do 3511 i＝1，igs，5

write（1，2993）file1（i），file1（i+1），file1（i+2），file1（i+3），file1（i+4）

2993 format（5（3x，a8，3x））

do 3778 iv＝1，ihv

f0＝fh（i，iv）

f1＝fh（i+1，iv）

f2＝fh（i+2，iv）

f3＝fh（i+3，iv）

f4＝fh（i+4，iv）

s0＝s（i，iv）

s1＝s（i+1，iv）

s2＝s（i+2，iv）

s3＝s（i+3，iv）

s4＝s（i+4，iv）

write（1，3556）iv，f0，s0，f1，s1，f2，s2，f3，s3，f4，s4

3778 continue

3556 format（i3，10f7.2）

3511 continue

close（1）

c 输出拟合节点计算的历时水位fh（i，iuv）、实测的历时水位sh（i，iuv）、拟合误差s（i，iuv），拟合节点名称file1（i）do 2013 i＝1，igs

tt＝0

filee＝file2（i）

open（1，file＝filee）

do 2018 iuv＝1，ihv-1

tt＝tt+tl（iuv）

write（1，2020）tt，fh（i，iuv），sh（i，iuv），s（i，iuv）

2018 continue

write（1，2030）tt，fh（i，iuv），sh（i，iuv），s（i，iuv），file1（i）

close（1）

2013 continue

2020 format（1x，f10.3，3f15.3，2x）

2030 format（1x，f10.3，3f15.3，2x，’＂’，a8，’＂’）

stop

end

8. 2. 1 模型的校验与识别

利用 1990～2000 年期间的水资源的实际利用量、河流水文、地下水位动态、气象等数据，对中游水资源数值模型进行校验与识别。

该期间实际水资源数据带入模型，模拟出地下水位动态过程、泉水流量过程、正义峡流量过程等模拟数据，将模拟数据与实际数据进行拟合对比，调整模型结构与模型参数，直至达到较好的拟合，即实现了对中游水资源系统的宏观模拟。

地下水位动态数据受地表水文随机因素、开采与灌溉随机性因素的影响，带有一定的随机性成分，某月的地下水位升降，或某季度甚至某年的水位变化趋势，并不一定能够代表区域地下水位的总趋势。因此某季度或某年的地下水下降值一般不能作为模型的校验依据。加之在数据处理中对实际水资源数据进行了一定的简化，简化归纳后的数据也带有微小的随机波动因素。

较理想的用于判断校对模型的数据，最好具有长时间系列、大变幅的特征，若变化幅度远超过随机波动干扰，其实质是将随机干扰 “过滤”掉了，提高模拟识别的可靠性。

中游干流平原区地下水研究程度较高，积累了大量的水利数据，选用 20 世纪 90 年代 10 年的地下水位累积变幅值与正义峡水文站历年 12 月至 2 月径流量数据作为模型校验识别依据 ( 12 月至2 月期间，由于黑河沿程不引水，此期间正义峡径流量基本上是泉水溢出量) ，尽管某年的地下水位动态具有一定的随机性，但累积 10 年的动态数据是非常可靠的。

经过调整模型结构参数与地层参数，使模拟水位与实际水位降深达到了较好的区域拟合，尤其是 10 年的累积水位变化量，相对拟合精度接近 90%。由此可以说明两方面的问题: 一是水资源数值模型的概化比较合理，二是该模型可较好地模拟水资源各要素之间的相互影响。经识别后的数值模型，可用于水资源调控预测及模拟分析，以科学合理配置中游地区的水资源。

8. 2. 2 地下水均衡分析

通过数值模拟，得到黑河干流中游平原区不同时期的地下水资源均衡结果 ( 表 8. 1) ，以及累计 10 年 ( 1990～1999 年) 地下水位降深模拟图 ( 图 8. 3) 。

图 8. 3 黑河干流中游平原区累计 10 年 ( 1990～1999 年) 地下水位降深模拟图

表 8. 1 中游平原地区地下水模拟均衡分析表 单位: 108m3/ a

由表 8. 1 分析可知，20 世纪 90 年代初、2000 年 ( 现状年) 两个典型时期的地下水总补给资源量分别为16. 627×108m3/ a 和 14. 632×108m3/ a，10 年间减少近 2×108m3/ a，其中渠系渗漏与田间灌溉渗漏减少 2. 125×108m3/ a。按整个中游平原计算区进行粗略统计，20 世纪 90 年代初干、支、斗渠的渠系利用系数约为 60%，田间灌溉入渗系数约为 15%～20%，到 2000 年，由于加强渠道防渗，干、支、斗三级渠系利用系数平均提高到 80%左右，地下水补给量大幅度减少，从而使地下水总补给量明显减少。

90 年代初至 2000 年这 10 余年间，为解决春旱问题，对地下水开采量有较大幅度提高，由0. 65×108m3/ a 逐步提高到 2. 476×108m3/ a，从而引起各地下水排泄要素重新调整，河水与泉水的溢出量及地下水蒸发量相应变小，河泉水溢出量由 90 年代初的 12. 131×108m3/ a 逐渐 减少 到11. 537×108m3/ a，减少了 0. 594×108m3/ a; 地下水蒸发量由 90 年代初的 4. 417×108m3/ a 逐渐减少到 2. 849×108m3/ a，减少了 1. 568×108m3/ a。

由图 8. 3 累计 10 年降深分布表明，地下水位降深大的位置，并没有大强度的地下水开采，显然不是开采地下水引起的。降深大的区域可超过 4m，最大值发生在民乐县洪水河与童子坝河山前的洪积扇上部，降深值超过十余米，其他降深大的位置，均沿南部山前埋深大且没有地下水开采的部位分布 ( 骆驼城地下水开采灌区除外) 。

在模型校正过程中，为寻求区域地下水降深的影响机制，对多种可能机制进行了大量组合模拟分析，经综合分析后得出结论: 产生如此形状降深场有两个主要的原因，其中最主要的原因是各灌区 “面状分布”渗漏量或灌溉回归补给量减少，即近十年来加强渠道衬砌防渗及逐步推广较省水灌溉方式形成的; 另一主要原因是山区拦蓄洪水使地下水山前补给量不断减少。

模拟结果同时表明，山前拦蓄洪水对地下水产生的后续影响将持续数十年甚至上百年才能达到新的平衡。

黑河是中游平原区最低的排泄基准面，在该种特定条件下，相对于泉水和蒸发排泄来说，河流溢出排泄量是相对稳定的，即增加地下水开采量，或者由于水利工程措施使地下水补给量减少，最先受到影响的应该是泉水上游的源头区溢出量与浅埋带地下水的蒸发量。由此，河流溢出量的衰减具有明显的滞后性，响应滞后周期长，而位于相对上游的泉水及浅埋带地下水蒸发，响应滞后周期较短，即泉水流量衰减相对较快。多年来的实际数据与模型模拟结论都证明了这一点，这与地下潜水的循环规律是相一致的。河流溢出量的大小，主要取决于河流附近的局部水力坡度，只要黑河附近地下水流场 ( 或坡度) 没有大的变化，河水溢出量就不会大幅度减少。当地下水埋深较小时，由于蒸发与埋深之间的非线性关系，地下水蒸发强度随地下水埋深急剧变化，虽然近十年来浅埋带地下水位下降幅度并不大，但地下水蒸发量却有较明显变化，尤其是在埋深小于 2m地区更为明显。当地下水位埋深超过 3m 后，降低地下水位所能夺取的地下水蒸发量有限。

从资源均衡的角度纵观中游干流平原地区地下水均衡，虽然整个计算区是负均衡的，但负均衡主要发生在远离黑河、泉及蒸发浅埋区的近山地带，具体表现为山前平原区地下水位的下降较多，黑河、泉及蒸发浅埋带水位降深小。以 2000 年均衡为例，在东南部 ( 民乐县) ，因地下水位持续下降而逐渐疏干上游区含水层，使该局部区域地下水负均衡量接近 1. 5×108m3/ a; 而靠近河流与泉水溢出带地区及地下水浅埋蒸发带，由于地下水排泄的 “自适应”调节作用 ( 当补给量减少时，排泄量将会自动缩减) ，地下水负均衡量较小，即在排泄带局部范围内，地下水补排大致平衡。

河流与泉水溢出量的响应滞后特征，容易给人们一种错觉，当某些水利工程运转之后，增加了部分地表水资源利用量，同时地下水补给量也随之减少。由于河流与或泉水响应滞后特性，其溢出量没有马上减少，表面上可利用的总水资源量 ( 地表水+地下水) 似乎增加了。这仅仅是短期的表现，实际情况是含水层 “地下水库”逐渐消耗，在较长的时期后，地下水溢出量减少会逐步表现出来，严重者使地下水资源枯竭。

以黑河中游平原东南部 ( 民乐县) 为例，地下水位比张掖附近的黑河水位高出 200 多米，当灌区地下水位下降不太大时，如 10m，相对于整个地下水位落差来说，其总体水力坡度变化还不到 10%，即在短期内，上游地带通过含水层向下游输送的地下水量不会明显减少 ( 短期内几乎是一个 “常数”) ，但要以不断疏干上游含水层为代价，据模拟均衡计算结果，现状条件下，每年疏干消耗民乐地区含水层地下水量约 ( 1. 5～2) ×108m3。从可持续发展的观点来看，长时期的疏干消耗上游含水层，一方面生态环境的极大改变不允许，同时将会导致地下水资源枯竭。这种开发利用方式可谓 “寅吃卯粮”，不能长时期持续。

随着渠道防渗工程的完善及节水技术的推广，使可利用的水量有所增加，应利用丰水年或丰水季节 “多余的”水资源对上游区进行回补，以阻止或减缓地下水资源向枯竭的方向演化，而不要盲目地扩大耕地面积，使水资源循环向合理可持续的状态转化。