Ntou123 -
3.2.2.1 水文地质概念模型
这是本节介绍的两种情况的第二种，即承压含水层既考虑其上覆定水头含水层的越流问题，又考虑其下伏定水头含水层的越流问题。上覆含水层、下伏含水层的水头简化为定水头，即整个计算期内不变，也不因中层水的开采而变。
与3.2.1节采用同一个水源地为例，其水文地质概念模型请参考3.2.1 节的概念模型，但计算时对条件的概化有所不同。
基本水文地质条件概化同于3.2.1节之处的是：非均质各向同性承压含水层、三边隔水、一边已知水头。
不同于3.2.1节之处的是：上覆的定水头潜水含水层通过弱透水层与计算的承压含水层有越流产生，下伏的定水头深层承压含水层通过弱透水层与计算的目的层也有越流产生，开采出的水量由四部分组成，第一部分是来自已知水头边界的侧向补给，第二部分是承压含水层的弹性释水量，第三部分是来自下伏的深层承压含水层的越流补给，第四部分是来自上部潜水含水层的越流补给。
据此建立二维流数学模型。
3.2.2.2 地下水流模型的建立
根据以上概化的水文地质概念模型，在考虑上覆潜水含水层、下伏深层承压含水层通过弱透水层对计算目的层越流补给的情况下，本区中层承压水可用下列非均质二维流数学模型来描述：
供水水文地质计算
式中：H1为H1（x，y）潜水含水层渗流区内各点任意时刻的水位（m）；H3为H3（x，y）深层承压含水层渗流区内各点任意时刻的水位（m）；T为中层承压含水层的导水系数（m2/d）；Qi为中层承压含水层中第i口井的开采量（m3/d）；δi为δ（x-xi，y-yi），为中层承压含水层中第i口井的δ函数，（xi，yi）为第i口井的坐标；S为中层承压含水层的弹性释水系数；H2（x，y，t）为中层承压含水层的水头函数（m）；H2，0（x，y）为中层承压含水层的初始水头函数（m）；B1，B2，B3为中层承压含水层的西、北、东边界，隔水边界；H2，1（x，y，x）为中层承压含水层南（B4）边界上的已知水头函数（m）；m为中层承压含水层抽水井的总数；C1为上层越流系数，C1＝K1′/M1′，K1′为上层弱透水层的渗透系数（m/d）；M1′为上层弱透水层的厚度（m）；C3为深层越流系数，C3＝K3′/M3′，K3′为深层弱透水层的渗透系数（m/d）；M3′为深层弱透水层的厚度（m）；D为计算区（82.68km2）。
以上数学模型可用有限元法求解。同3.2.1将计算区剖分成108个三角形单元，见图3-2-1；剖分成67个节点，其中未知水头节点59个，已知水头的一类边界节点8个，未知水头的二类边界节点16个，边界节点共24个，内节点43个，见图3-2-2。
3.2.2.3 计算程序
把求解以上数学模型的Fortran有限单元源程序复制如下。该程序是用来调参的，所以其时段个数的配置、数据文件的数组配置是与调参所用的原始数据完全匹配的。该程序是在3.2.1程序的基础上又增添了潜水、深层承压含水层越流有关内容后形成的，与3.2.1的程序大部分语句是相同的。
行号 源 程 序
1 c f3.for 承压水有上（潜）、下（深承）两层越流（对弱透水层只考虑越流不考虑弹性释放），调参用
2 c 上（潜）、下（深承）两层不剖分，借用中层的节点号只有定水位，靠水头差产生越流
3 c 6 行：可调数组的各个参数预先赋值，不同计算区只要改变这些参数即可用此程序，其中：
4 c ms：单元个数，ids：节点个数；n：未知水头节点数；mi：开采节点数；igs：拟合点数
5 c icqs：参数区个数，ihv：计算时段数，nn：已知水头节点个数
6 parameter（ms＝108，ids＝67，n＝59，mi＝30，igs＝14，icqs＝8，ihv＝24，nn＝8）
7 c h0：时段初水头（m）；hed：时段末水头（m）；in：各单元三节点编号（必须按小中大顺序排），fh：拟合点各时段计算水位（m）；
8 dimension h0（ids），h01（ids），h03（ids），hed（ids），in（ms，3），fh（igs，ihv）
9 c sh：拟合点各时段实测水位（m）；idyh：导水矩阵工作单元 zb：各节点 x，y坐标（输入坐标为图面mm数）；igdh：拟合点节点号；q：开采井各时段水量（抽为正，注为负，流出边界为正，流入为负）（m3/d）
10 dimension sh（igs，ihv），idyh（ids，9），zb（ids，2），igdh（igs），q（mi，ihv）
11 c d：各节点导水矩阵；ifdh：开采井节点号；s：拟合点拟合误差（m）；mqh：各单元的参数区号；e：各节点释（储）水矩阵；
12 dimension ifdh（mi），d（ids，9），s（igs，ihv），mqh（ms），e（ids）
13 c fqc：各参数区最多4个参数之值；sss：各节点水头总变幅，由 t0时刻起算（m）；h：初始流场，yo1：各节点潜越流矩阵工作单元，yo3：各节点向深层越流矩阵工作单元，tl：各时段步长值（天）；bc：三角单元几何量（bi，ci，bj，cj，bk，ck）；h0nn：已知变水头节点各时段的水头
14 dimension fqc（icqs，4），sss（ids），h（ids），yo1（ids），yo3（ids），tl（ihv）
15 dimension bc（3，2），h0nn（nn，ihv）
16 c file1（igs）*8，file2（igs）*8：定义8字符的字符串各 igs个；filee*8：定义8字符的文件名
17 character file1（igs）*8，file2（igs）*8，filee*8
18 write（*，23）
19 23 format（1x，"zz1＝？ zz2＝？ 迭代因子：zz1潜水0—1，zz2承压水1—2"）
20 c zz2：承压水为超松弛系数，当 ids＜300时取1.2～1.3，当 ids＞300时取1.3～1.5为好
21 zz2＝1.25
22 write（*，24）ms，ids，n，mi，igs，icqs，ihv
23 24 format（5x，12i5）
24 open（1，file＝"fqc"，status＝"old"）
25 c从"fqc"中读取承压水的区号m10（空读），参数1（T），参数2（μ*），参数3潜-承的越流系数（k"/m"），参数4：中-深的越流系数（k"/m"）
26 read（1，*）（m10，fqc（i，1），fqc（i，2），fqc（i，3），fqc（i，4），i＝1，icqs）
27 close（1）
28 c"in"：各单元三节点编号文件（必须按小中大顺序排），ia为单元号空读
29 open（1，file＝"in"，status＝"old"）
30 read（1，*）（ia，（in（i，j），j＝1，3），i＝1，ms）
31 close（1）
32 c igdh：拟合点节点编号
33 open（1，file＝"igdh"，status＝"old"）
34 read（1，*）（igdh（i），i＝1，igs）
35 close（1）
36cmqh：各单元的参数区号，ia为单元号空读
37 open（1，file＝"mqh"，status＝"old"）
38 read（1，*）（ia，mqh（i），i＝1，ms）
39 close（1）
40 c ifdh：开采井节点号
41 open（1，file＝"ifdh"，status＝"old"）
42 read（1，*）（ifdh（i），i＝1，mi）
43 close（1）
44 c tl（i）：各时段值（天/时段，每月为1个时段）
45 open（1，file＝"tl1"，status＝"old"）
46 read（1，*）（tl（i），i＝1，ihv）
47 close（1）
48 c"zb"：各节点 x，y坐标文件（此处为1:2.5万图面的mm数），ia为节点号空读
49 open（1，file＝"zb"，status＝"old"）
50 read（1，*）（ia，（zb（i，j），j＝1，2），i＝1，ids）
51 close（1）
52 c把各节点 x，y坐标1:2.5万图面的mm数换算为实地的m数
53 do 888 i＝1，ids
54 do 888 j＝1，2
55 zb（i，j）＝zb（i，j）*25
56 888 continue
57 c中层水开采井第1时段的开采量（m3/d），空读 ia开采井个数，空读 ib开采井所在节点编号
58 open（1，file＝"qz1"，status＝"old"）
59 read（1，*）ia，（ib，q（i，1），i＝1，mi）
60 close（1）
61 c中层水开采井第13时段的开采量（m3/d），空读 ia开采井个数，空读 ib开采井所在节点编号
62 open（1，file＝"qz2"，status＝"old"）
63 read（1，*）ia，（ib，q（i，13），i＝1，mi）
64 close（1）
65 c中层水开采井第2～12时段的开采量等于第1时段的开采量
66 do 887 j＝2，12
67 do 886 i＝1，mi
68 q（i，j）＝q（i，1）
69 886 continue
70 887 continue
71 c中层水开采井第14-ihv时段的开采量等于第13时段的开采量
72 do 884 j＝14，ihv
73 do 883 i＝1，mi
74 q（i，j）＝q（i，13）
75 883 continue
76 884 continue
77 c sh：拟合点各时段实测水位（m），空读 ia拟合点所在节点编号
78 open（3，file＝"sh"，status＝"old"）
79 read（3，*）（ia，（sh（i，j），j＝1，ihv），i＝1，igs）
80 close（3）
81 c qc：潜水初始流场-定水位流场，空读 is节点编号
82 open（1，file＝"qc"，status＝"old"）
83 read（1，*）（is，h01（i），i＝1，ids）
84 close（1）
85 c zc：中层水初始流场，空读 is节点编号
86 open（1，file＝"zc"，status＝"old"）
87 read（1，*）（is，h（i），i＝1，ids）
88 close（1）
89 c sc：深层水初始流场-定水位流场，空读 is节点编号
90 open（1，file＝"sc"，status＝"old"）
91 read（1，*）（is，h03（i），i＝1，ids）
92 close（1）
93 c h0nn：中（h0nn）层已知变水头节点各时段的水头，空读 ii节点编号
94 open（1，file＝"h0nn"）
95 read（1，*）（ii，（h0nn（i，ikv），ikv＝1，ihv），i＝1，nn）
96 close（1）
97 c"file1"存放 igs个字符串（拟合节点名称占4字符，及观测井的原编号再占4字符）
98 open（1，file＝"file1"，status＝"old"）
99 read（1，110）（file1（i），i＝1，igs）
100 close（1）
101 c"file2"存放 igs个字符串（拟合节点名称占4字符，加.dat后缀再占4字符）
102 open（1，file＝"file2"，status＝"old"）
103 read（1，110）（file2（i），i＝1，igs）
104 close（1）
105 110 format（10a8）
106 c ：计算开始前先把全部节点的时段末刻水头 hed（i），时段初刻水头 h0（i）赋初值 h（i），以使开始迭代时 hed（i），h0（i）不为零
107 do 1993 i＝1，ids
108 hed（i）＝h（i）
109 1993 h0（i）＝h（i）
110 c对中层承压水导水矩阵工作单元 idyh，释水矩阵 e，越流矩阵 yo1，yo3，导水矩阵 d，赋0；
111 do 25 i＝1，ids
112 do 25 j＝1，9
113 idyh（i，j）＝0
114 e（i）＝0.0
115 yo1（i）＝0.0
116 yo3（i）＝0.0
117 25 d（i，j）＝0.0
118 c sum2用来累计计算区总面积，先赋零
119 sum2＝0
120 c对承压水逐个单元计算几何量及导水、释水、越流矩阵；对承压水第 ip单元的三个节点号依次赋给 i，j，k及i1，j1，k1
121 do 80 ip＝1，ms
122 i＝in（ip，1）
123 j＝in（ip，2）
124 k＝in（ip，3）
125 i1＝i
126 j1＝j
127 k1＝k
128 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中的位置
129 c且先占者排前，94循环体：使 j1，k1分两轮到 idyh中找位置；当 j1找时，134行发现 i1的 idyh中已有 j1，则跳到 j1对 i1的94循环体头，换为 k1对 i1循环；当134行没发现已有 j1，且135行判断此位不空时返回96循环体头找下一个位置，当136行碰到第1个空位时，由 j1占据，然后跳到 j1对 i1循环体96头，换为 k1对 i1循环
130 do 90 iv＝1，3
131 idyh（i1，1）＝i1
132 do 94 iu＝j1，k1，k1-j1
133 do 96 i2＝2，9
134 if（idyh（i1，i2）.eq.iu）goto 94
135 if（idyh（i1，i2）.ne.0）goto 96
136 if（idyh（i1，i2）.eq.0）idyh（i1，i2）＝iu
137 goto 94
138 96 continue
139 94 continue
140 c把 ip单元的下一个节点 j提为 i1，k提为 j1，i降为 k1，然后返回90循环头处理 ip单元的下一个节点，循环3次则 ip单元中 i，j，k 3个点的 idyh全部找完
141 iuu＝i1
142 i1＝j1
143 j1＝k1
144 k1＝iuu
145 90 continue
146 c第 ip单元中 i，j，k三节点的 x坐标赋给 xi，xj，xk，y坐标赋给 yi，yj，yk
147 xi＝zb（i，1）
148 xj＝zb（j，1）
149 xk＝zb（k，1）
150 yi＝zb（i，2）
151 yj＝zb（j，2）
152 yk＝zb（k，2）
153 c第 ip单元三角形面积 ss
154 ss＝abs（（xj-xi）*（yk-yi）-（xk-xi）*（yj-yi））*0.5
155 c累加各单元面积
156 sum2＝sum2+ss
157 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）
158 bc（1，1）＝yj-yk
159 bc（1，2）＝xj-xk
160 bc（2，1）＝yk-yi
161 bc（2，2）＝xk-xi
162 bc（3，1）＝yi-yj
163 bc（3，2）＝xi-xj
164 c第 ip单元所在参数区号赋给 jv，参数1（T）赋给 txy，参数2（μ*）赋给 ts，参数3（k"/m"）赋给 rmk1，参数4（k"/m"）赋给 rmk3
165 jv＝mqh（ip）
166 txy＝fqc（jv，1）
167 ts＝fqc（jv，2）
168 rmk1＝fqc（jv，3）
169 rmk3＝fqc（jv，4）
170 c第 ip单元三节点 i，j，k的释水矩阵元素 c（ii，ii）及 c（ii，jj）**式中没包括1/（2Δt）时，随着 ip＝1～ms的单元循环而对有关单元求其和
171 e（i）＝-ts*ss/3.0+e（i）
172 e（j）＝-ts*ss/3.0+e（j）
173 e（k）＝-ts*ss/3.0+e（k）
174 c越流矩阵元素，第 ip单元1m水头差时的越流量（m3/d/m）均分给三节点 i，j，k，随着 ip＝1～ms的单元循环而对有关单元求其和
175 yo1（i）＝yo1（i）+rmk1*ss/3.0
176 yo1（j）＝yo1（j）+rmk1*ss/3.0
177 yo1（k）＝yo1（k）+rmk1*ss/3.0
178 yo3（i）＝yo3（i）+rmk3*ss/3.0
179 yo3（j）＝yo3（j）+rmk3*ss/3.0
180 yo3（k）＝yo3（k）+rmk3*ss/3.0
181 c计算 ip单元各节点的导水矩阵元素 d（i，j），（i＝1，2，3）（j＝1，2，3）
182 do 100 iu＝1，3
183 do 104 iuu＝1，3
184 c当 iu＝1时，即 ip单元的 i节点，分别计算 iuu＝1，2，3，即 ip单元的 i，j，k节点对 i节点的 ai，ai即公式**的没求和部分
185 i＝in（ip，iu）
186 j＝in（ip，iuu）
187 ai＝txy*（bc（iu，1）*bc（iuu，1）+bc（iu，2）*bc（iuu，2））/ss/4.0
188 c在 ip单元的三个节点中，排除第3点，只让第1，2点（i，j点）的 ai加入 i点的导水矩阵元素 d（i，j）中，第186行的 j可分别轮到 i，j，k三点，但第189行的1～9个中，仅有 i，j点在 i点的 idyh中，此句排除了第3点加入 i点的导水矩阵元素 d（i，j）中随着 ip＝1-ms的循环，到194行时承压水导水矩阵已完全形成
189 do 106 k＝1，9
190 if（j.eq.idyh（i，k））d（i，k）＝ai+d（i，k）
191 106 continue
192 104 continue
193 100 continue
194 80 continue
195 c至此，中层几何量计算完，以下开始时段循环
196 c在时段循环中，承压水与时段有关的部分在时段循环中完成
197 do 280 ikv＝1，ihv
198 write（*，*）"ikv＝"，ikv
199 c该时段步长赋给 dt，最大误差记录 amax，开采井点录用计数器 iqq
200 dt＝tl（ikv）
201 iqq＝1
202 988 amax＝0.0
203 do 232 i＝1，nn
204 hed（n+i）＝h0nn（i，ikv）
205 h0（n+i）＝h0nn（i，ikv）
206 232 continue
207 do 400 i＝1，n
208 c其中的潜水越流流入量采用（本节点潜水定水位 h01（i）减承压水本时段末水位 hed（i））之差 hedd计算
209 hedd＝h01（i）-hed（i）
210 c其中的深层水越流流出量采用（本节点承压水本时段末水位 hed（i）减本节点深层水定水位 h03（i））之差heddd计算
211 heddd＝hed（i）-h03（i）
212 res＝-hedd*yo1（i）+heddd*yo3（i）-e（i）*（hed（i）-h0（i））/dt
213 50 if（i.eq.ifdh（iqq））then
214 res＝res+q（iqq，ikv）
215 iqq＝iqq+1
216 endif
217 do 300 j＝1，9
218 k＝idyh（i，j）
219 if（k.eq.0）goto 300
220 res＝res+d（i，j）*hed（k）
221 300 continue
222 res＝res/（-d（i，1）+e（i）/dt）
223 hed（i）＝res*zz2+hed（i）
224 if（abs（res）.le.amax）goto 400
225 amax＝abs（res）
226 imax2＝i
227 400 continue
228 write（*，*）"amax2＝"，amax，imax2
229 if（amax.gt.9.999999e-03）goto 988
230 c至此，中层水全部节点末刻水头迭代完一轮，节点误差最大者如果 ＞0.01m，则返回 988句进行下一轮迭代，如果≤0.01m，则迭代结束，中层水本时段末刻流场计算完
231 c以下开始把本时段末刻水位中层水全部节点本时段末刻水头 hed赋给下时段初刻水位 h0，拟合点本时段末刻水位 hed存入 fh
232 189 do 281 i＝1，ids
233 h0（i）＝hed（i）
234 281 continue
235 do 190 j＝1，igs
236mp＝igdh（j）
237 fh（j，ikv）＝hed（mp）
238 190 continue
239 280 continue
240 c至此，时段循环完
241 c屏幕输出计算区总面积 sum2（m2）
242 write（*，*）"sum＝"，sum2，"m2"
243 do 1926 i＝1，ids
244 1926 sss（i）＝hed（i）-h（i）
245 c"ok"输出中层水的最终流场：节点号 i、坐标 zb、最终流场 hed、该节点0-ihv时间（最终流场比初始流场）的水位总变幅
246 open（1，file＝"ok"）
247 write（1，1927）（i，zb（i，1），zb（i，2），hed（i），sss（i），i＝1，ids）
248 close（1）
249 1927 format（2x，i4，2f9.0，2f9.2）
250 c计算出拟合点各时段误差（计算水位-实测水位）s（i，iuv）
251 do 8 iuv＝1，ihv
252 do 2003 i＝1，igs
253 s（i，iuv）＝fh（i，iuv）-sh（i，iuv）
254 2003 continue
255 8 continue
256 c输出拟合点各时段误差
257 open（1，file＝"gan"）
258 do 3511 i＝1，igs，7
259 write（1，*）
260 write（1，2993）file1（i），file1（i+1），file1（i+2），file1（i+3），file1（i+4），
261 * file1（i+5），file1（i+6）
262 2993 format（1x，7（3x，a8，3x））
263 do 3778 iv＝1，ihv
264 f0＝fh（i，iv）
265 f1＝fh（i+1，iv）
266 f2＝fh（i+2，iv）
267 f3＝fh（i+3，iv）
268 f4＝fh（i+4，iv）
269 f5＝fh（i+5，iv）
270 f6＝fh（i+6，iv）
271 s0＝s（i，iv）
272 s1＝s（i+1，iv）
273 s2＝s（i+2，iv）
274 s3＝s（i+3，iv）
275 s4＝s（i+4，iv）
276 s5＝s（i+5，iv）
277 s6＝s（i+6，iv）
278 write（1，3556）iv，f0，s0，f1，s1，f2，s2，f3，s3，f4，s4，f5，s5，f6，s6
279 3778 continue
280 3556 format（i2，1x，7（f7.2，f5.2，2x））
281 3511 continue
282 close（1）
283 c输出拟合节点计算的历时水位 fh（i，iuv）、实测的历时水位 sh（i，iuv）、拟合误差 s（i，iuv），拟合节点名称 file1（i）
284 do 2013 i＝1，igs
285 tt＝0
286 filee＝file2（i）
287 open（1，file＝filee）
288 do 2018 iuv＝1，ihv-1
289 tt＝tt+tl（iuv）
290 write（1，2020）tt，fh（i，iuv），sh（i，iuv），s（i，iuv）
291 2018 continue
292 write（1，2030）tt，fh（i，iuv），sh（i，iuv），s（i，iuv），file1（i）
293 close（1）
294 2013 continue
295 2020 format（1x，f10.3，3f15.3，2x）
296 2030 format（1x，f10.3，3f15.3，2x，"″"，a8，"″"）
297 stop
298 end
3.2.2.4 模型校正所需数据
这里所需数据在3.2.1.4中已大部分列出，除此之外需增加的数据：一是分区参数文件fqc需再增补一列，即每区再增补一个越流系数；增加的数据二是潜水含水层的定水位流场qc；增加的数据三是深层承压水含水层的定水位流场sc。下面分别copy 列出。
一、fqc（注：分栏排列）
1 1000 0.008 0.00001 0.00001 5 1700 0.0005 0.0001 0.0001
2 3500 0.0001 0.0005 0.0005 6 3500 0.0005 0.000001 0.000001
3 3500 0.0005 0.0007 0.0007 7 1900 0.0005 0.00001 0.00001
4 3500 0.0005 0.0007 0.0007 8 2000 0.0005 0.00001 0.00001
二、qc
01 379.50 02 378.70 03 378.20 04 378.50 05 376.87
06 378.50 07 375.90 08 374.55 09 378.50 10 377.50
11 378.50 12 377.77 13 378.50 14 379.50 15 380.10
16 379.00 17 379.50 18 379.00 19 381.00 20 379.10
21 378.90 22 378.80 23 377.70 24 377.90 25 377.10
26 376.80 27 375.30 28 375.38 29 374.56 30 374.90
31 375.20 32 379.70 33 381.20 34 383.20 35 381.80
36 381.33 37 380.00 38 375.60 39 377.40 40 376.70
41 377.10 42 376.20 43 377.77 44 378.58 45 378.82
46 380.00 47 382.50 48 382.00 49 379.90 50 380.60
51 380.83 52 378.50 53 382.60 54 377.80 55 378.60
56 382.09 57 382.20 58 382.10 59 382.20 60 383.00
61 382.60 62 380.30 63 379.60 64 380.20 65 382.00
66 382.86 67 381.10
三、sc
1 372.6 2 372.6 3 371.5 4 372.6 5 371.5
6 372.6 7 371.4 8 371.3 9 372.2 10 370.5
11 372.1 12 370.73 13 372.1 14 372.1 15 370.5
16 370