大区域地下水模拟的理论和方法[J].doc

22大区域地下水模拟的理论和方法张祥伟（中国水利水电科学研究院 北京 100044）【摘 要】 大区域地下水数值模拟中，存在时空尺度选择和信息不足问题本文在融合地质统计、逆问题理论和地下水运动理论的基础上，提出了建立大区域地下水运动模拟的理论和方法主要内容包括：（1）提出了大区域地下水有限元计算的时空尺度理论公式；（2）针对信息不足，提出了大区域地下水位推定的方法；（3）根据推定的地下水位，运用逆解析理论法对大区域地下水含水层透水系数进行逆推定；（4）运用实例关键词】 大区域地下水，数值计算，逆问题1. 引言在我国许多地区，由于水资源严重短缺，超采地下水，导致大范围的地下水位急剧下降，从而引起区域性的地面沉降、海水入侵以及湿地退化等生态环境问题，严重影响了区域经济社会可持续发展如何合理地利用地下水资源，对区域和流域地下水、地面水进行优化配置，通过建立健全的流域水循环体系，保护地下水环境，恢复流域生态，实现水资源的可持续利用以及人与自然的协调，是当今水资源研究领域十分重要的课题要合理地利用区域地下水资源，首要问题是要对区域地下水资源进行符合实际的模拟，摸清区域地下水的运动规律然而，由于区域地下水模拟涉及大尺度地下水运动，而传统的地下水模拟手段和方法在处理大尺度地下水运动模拟中所涉及的有限元计算中的时空尺度的选择、水位和含水层信息不足等问题上方面存在困难，引起由于信息不足带来的较大计算误差。

因此，针对尺度和信息不足问题，建立大区域地下水运动模拟的理论和方法十分必要本文在融合地质统计、逆问题理论和地下水运动理论的基础上，提出了建立大区域地下水运动模拟的理论和方法首先，提出了大区域地下水有限元计算的时空尺度理论公式；其次，针对信息不足，提出了大区域地下水位推定的方法；第三，根据推定的地下水位，运用逆解析理论法对大区域地下水含水层透水系数进行逆推定；最后，给出了运用实例2. 大区域地下水有限元计算的时空尺度在大区域地下水数值计算中，首先需要选择研究对象所需划分的计算网格大小或确定计算时间步长大小尤其在模拟河道网基础上的分布式水文模型和大区域地下水模型进行地表水与地下水偶合计算时，同时兼顾河道网的数值地形网格大小和地下水网格大小，既满足水文计算的精度同时使地下水计算不至于失稳就显得至关重要然而，由于大区域地下水计算23中，由于范围较大，人们往往希望选择较大的网格进行计算，以减少计算工作量这往往导致计算上的不稳定对于地下水有限元计算的时空尺度（时空步长）的选择，1977 年，Newman 等对二维地下水运动方程的有限元解法中的不稳定问题进行了分析 [1]，推测不合理的时间步长导致了混合型差分的出现，由此导致了解的不稳定问题。

针对 Neman 等的推测，Zhang （1992） [2] 和Wood（1996） [3]提出了二维地下水有限元计算的时间步长的条件2002 年，张祥伟等 [4]运用最大最小原理对大尺度二维和准三维地下水有限元计算的不稳定问题进行了理论分析，提出了二维和准三维地下水有限元计算的时空步长的理论公式即对于准三维地下水有限元计算的时空步长分别为：空间步长条件：（1）),min(2188axmSktTktL其中，L max 为准三维地下水有限元计算的三角形网格中的最小三角形最小边长； Δt 为时间步长； 、 分别为层间难透水层的厚度和透水系数；T 1、T 2 为非承和承压含水层的透水量mk系数；μ、S 分别为非承和承压含水层的储留系数二维的情况下，空间步长的条件为：（2）StTL83. 大区域地下水位的推定在进行大区域地下水计算中，需要根据实测的地下水位推定初始流场，以便检验地下水数值计算精度、进行非恒定地下水计算以及识别含水层参数对于初始流场的推定的方法，通常有 Universal Kriging 方法[5-6] ，UK 法的计算行列比较大，计算比较复杂Newman 等（1984） [7] 和 Sun（1999） [8] 运用 Residual Kriging（RK）法进行区域地下水位的推定，也就是将实测地下水趋势面去除得到正态残差，将正态残差运用 Ordinary Kriging 法进行面上残差的推定，再加上实测地下水面的趋势得到三角形网格上各点的推定地下水位值。

上述方法在大区域地下水计算中遇到信息不足的问题，当实测地下水位信息不足时，推定的初始流场会带来较大的误差，本文根据地形水文学的原理[9]，即地下水与地形之间存在的相关性，提出运用数值地形模型（DEM）中的地形标高作为辅助信息，修正实测地下水位得到的地下水趋势面，然后，运用 Ordinary Kriging 方法对修正后趋势面得到的正态残差对三角形网格点上的残差进行推定，每各网格点上推定的残差再加上修正的趋势值得到面上地24下水位的推定值本文将此方法定名为 ROKMT（Residual Ordinary Kriging with Modified Trend）方法 [10-11]，即地下水位由以下两部分构成：（3）()()hmxx其中， 为实测地下水位的趋势面残差， 为 DEM 修正后的地下水趋势面，即() ()mxyazya5243210)(（4）98276yaxz其中，a 0-a9 为趋势方程系数， x、y 为坐标，z 为 DEM 的地形标高当残差 满足正态性条件时，用 Ordinary Kriging 法进行推定，再加上式（4）即可()计算得到大区域地下水的推定值大区域地下水位推定的方法如图 1。

4. 区域地下水透水系数的识别在地下水流动模拟中，一般根据含水层的水文地质条件和土壤特性，选定透水系数的取值在大区域地下水计算中，含水层物理特性的信息常常十分有限，即使知道含水层的土壤类型，但由于同类土壤的透水系数的变化范围很大，最小值与最大值之间甚至有 100 倍以上的变化幅度这为正确选择地下水参数带来困难，本文针对式（5）的二维地下水运动，根据推定的地下水位值，运用 Guass-Newton 法对大区域地下水含水层透水系数进行识别二维地下水运动的基本方程为：（5）()()0hhTqxy其中，h 为地下水位，T 为透水量系数，q 为地下水涵养量或扬水量方程式（5）的透水量系数的识别，有直接法和间接法直接法是根据已知的地下水位值，Yes 、、 、OK、、、、、 、 、、 入力 、、 、 No 、 、 、、 + 、、、 、1 、 、、 、 、 、、(4)、、、、、Kolmgorv -Smirnov、、、、 25直接从式（5）中反求透水量系数 T直接法的问题是，由于地下水位中含有误差，细小的误差将导致水位的偏微分很大的误差，计算稳定性差，而且往往引起不适定问题，大大降低了参数的计算精度。

本文运用间接法识别水文地质参数即给定含水层透水量 T 的初始值进行迭代计算，根据 3 中推定的地下水位值和反复计算得到的计算水位值的残差平方和最小作为目标函数，计算最优的透水量值 T目标函数为：（6）1211(,)[()()]kmkkobsl Thf h其中，h obs为 ROKMT 法的推定地下水位值 Tk+1为第 K 次迭代的透水量系数通过下式计算，（7）1()kjk obsThTh根据含水层的物理特性，透水量系数需满足以下限制条件：（8）minaxjjj其中，T jminhe和 Tjmax为透水量系数的下限和上限值当透水量系数识别完成后，运用不规则三角形差分（TFDM）法进行大区域地下水计算综合 3 和 4 的计算步骤，大区域地下水计算方法的如图 25. 运用实例-Sarobetsu 湿地地下水模拟本文将所提出的方法运用到日本北海道 Sarobetsu 湿地地下水模拟中5.1 研究对象的概况、 、、、、 、、 、、、、、 、、、、 、 、、、、 、 、、、、、 、、 、、、 TFDM、 、、、、、、、 ROKMT、 、、、、 、DEM、、、 、、、、、、、、 、 、、 、 、 2 、、 、、、 26Sarobetsu 湿地位于日本最北端（图 3） ，为日本最大的以水台癣为主的高层湿地，面积635km2。

近年来，由于农田开垦等土地开发利用，地下水位降低，湿地出现退化，海水上朔，生态环境恶化最主要的表现是：5.1.1 湿地指示性植物减少对地下水位敏感的湿地植物水台癣面积逐渐减少，低竹类植物风长，面积逐年扩大，植物耗水量增加，植物叶面蒸散发增强5.1.2 湖泊面积减小由于湿地干旱化导致 1993 年湿地内三大湖泊兜湖、Paken 湖和Peken 湖的面积比 1975 年减少 10%三湖泊周边的水台癣急剧减少，生态环境恶化5.1.3 地面沉降湿地中心区地面最大沉降达到 50 厘米5.1.4 海水上朔由于地面沉降造成和 Sarobetsu 河的枯水期流量减少，海水沿Sarobetsu 河口上朔流入 Paken 湖和 Peken 湖，破坏了湖周边生态图 3 Sarobetsu 湿地位置图 图 4 Sarobetsu 地形图5.1.5 Sarobetsu 河下游洪水泛滥 由于 Sarobetsu 湿地的地面沉降和海水的顶托，造成下游洪水泛滥因此，为恢复湿地和减少地面沉降，首要课题是研究该区域地下水位恢复的措施前提是研究该区域地下水涵养机理、地下水的运动规律。

5.2 地形和地质Sarobetsu 湿地的地形为从北到南沿海岸线逐渐降低Sarobetsu 河和海岸之间为 5-20 米N研究区域日本海清明川兜沼川天盐川福永川日本北海道Sarobetsu 河Paken 湖Peken 湖兜 湖05015020Distance(m)501502503Distance (m)丘陵低地日本海ABB1C1CA1地质调查测线27的沙丘地带，湿地地面高程在 2 米至 8 米之间，周边为起伏的台地和丘陵（如图 4） 该区域地质上主要分为两层，即更别层以上为泥岩层，为由南向北的向斜构造，从东、北、西向日本海方向嵌入更别层以上为 30-40 米的洪积层，主要以透水性强的泥炭和沙质层为主更别层以下为化石地下水，与上部地下水交换很少本研究以更别层以上低地的非承压地下水运动为对象5.3 地下水位空间分布的推定运用本文提出的 ROKMT 法对 Sarobetsu 湿地地下水位进行推定5.3.1 网格的划分基于 250x250m 的 DEM 地形标高扩展成 500x500m 的矩形网格，每个矩形网格分为两格三角形全对象区域分为 1051 个节点，1903 个三角形网格（如图 5） 。

5.3.2 地下水位的推定结果根据 1997 年 72 个地下水监测点同步地下水监测值、5 个 Sarobetsu 河的水位值和 3 个湖泊的水位值，运用 ROKMT 法推定地下水位另外的 10 个实测地下水位作为验证点（图 6） 地下水位的推定结果如图 710 个验证点的平均误差（ME） ，平均绝对误差（MAE）和平方根二乘误差（RMSE）分别为 0.005m，0.437m 和 0.598m从三种误差看，推定的地下水位空间分布具有较好的精度，可以用来识别透水量系数5.3.3 透水系数的识别结果（1）边界条件的设定海岸线一侧和研究区内部的河流、湖泊作为定水头边界其他边界位于丘陵的边沿作为流量边界2）涵养率（入渗补给系数）和透水量系数的界限值的设定兜湖PakenPeken地下水位监测点河流湖泊124569 871030 50001～10: 验证点28图 5 计算网格 图 6 地下水位验证点位置图图 7 地下水位的推定结果 图 8 透水系数识别结果运用分布式 Tank 水文模型，对 1997 年降雨径流的模拟，以及水平衡法水量平衡计。