大区域地下水模拟的理论和方法

　1. 引言

　　在我国许多地区，由于水资源严重短缺，超采地下水，导致大范围的地下水位急剧下降，从而引起区域性的地面沉降、海水入侵以及湿地退化等生态环境问题，严重影响了区域经济社会可持续发展。如何合理地利用地下水资源，对区域和流域地下水、地面水进行优化配置，通过建立健全的流域水循环体系，保护地下水环境，恢复流域生态，实现水资源的可持续利用以及人与自然的协调，是当今水资源研究领域十分重要的课题。
　　要合理地利用区域地下水资源，首要问题是要对区域地下水资源进行符合实际的模拟，摸清区域地下水的运动规律。然而，由于区域地下水模拟涉及大尺度地下水运动，而传统的地下水模拟手段和方法在处理大尺度地下水运动模拟中所涉及的有限元计算中的时空尺度的选择、水位和含水层信息不足等问题上方面存在困难，引起由于信息不足带来的较大计算误差。因此，针对尺度和信息不足问题，建立大区域地下水运动模拟的理论和方法十分必要。
　　本文在融合地质统计、逆问题理论和地下水运动理论的基础上，提出了建立大区域地下水运动模拟的理论和方法。首先，提出了大区域地下水有限元计算的时空尺度理论公式；其次，针对信息不足，提出了大区域地下水位推定的方法；第三，根据推定的地下水位，运用逆解析理论法对大区域地下水含水层透水系数进行逆推定；最后，给出了运用实例。
　　2. 大区域地下水有限元计算的时空尺度
　　在大区域地下水数值计算中，首先需要选择研究对象所需划分的计算网格大小或确定计算时间步长大小。尤其在模拟河道网基础上的分布式水文模型和大区域地下水模型进行地表水与地下水偶合计算时，同时兼顾河道网的数值地形网格大小和地下水网格大小，既满足水文计算的精度同时使地下水计算不至于失稳就显得至关重要。然而，由于大区域地下水计算中，由于范围较大，人们往往希望选择较大的网格进行计算，以减少计算工作量。这往往导致计算上的不稳定。
　　对于地下水有限元计算的时空尺度（时空步长）的选择，1977年，Newman等对二维地下水运动方程的有限元解法中的不稳定问题进行了分析[1]，推测不合理的时间步长导致了混合型差分的出现，由此导致了解的不稳定问题。针对Neman等的推测，Zhang（1992）[2] 和Wood（1996）[3]提出了二维地下水有限元计算的时间步长的条件。2002年，张祥伟等[4]运用最大最小原理对大尺度二维和准三维地下水有限元计算的不稳定问题进行了理论分析，提出了二维和准三维地下水有限元计算的时空步长的理论公式。即对于准三维地下水有限元计算的时空步长分别为：
　　空间步长条件：
　　
　　3. 大区域地下水位的推定
　　在进行大区域地下水计算中，需要根据实测的地下水位推定初始流场，以便检验地下水数值计算精度、进行非恒定地下水计算以及识别含水层参数。
　　对于初始流场的推定的方法，通常有Universal Kriging方法[5-6]，UK法的计算行列比较大，计算比较复杂。Newman等（1984）[7] 和Sun（1999）[8] 运用Residual Kriging（RK）法进行区域地下水位的推定，也就是将实测地下水趋势面去除得到正态残差，将正态残差运用Ordinary Kriging法进行面上残差的推定，再加上实测地下水面的趋势得到三角形网格上各点的推定地下水位值。
　　上述方法在大区域地下水计算中遇到信息不足的问题，当实测地下水位信息不足时，推定的初始流场会带来较大的误差，本文根据地形水文学的原理[9]，即地下水与地形之间存在的相关性，提出运用数值地形模型（DEM）中的地形标高作为辅助信息，修正实测地下水位得到的地下水趋势面，然后，运用Ordinary Kriging方法对修正后趋势面得到的正态残差对三角形网格点上的残差进行推定，每各网格点上推定的残差再加上修正的趋势值得到面上地下水位的推定值。本文将此方法定名为ROKMT（Residual Ordinary Kriging with Modified Trend）方法[10-11]，即地下水位由以下两部分构成：
　　
　　大区域地下水位推定的方法如图1。
　　
　　4. 区域地下水透水系数的识别
　　在地下水流动模拟中，一般根据含水层的水文地质条件和土壤特性，选定透水系数的取值。在大区域地下水计算中，含水层物理特性的信息常常十分有限，即使知道含水层的土壤类型，但由于同类土壤的透水系数的变化范围很大，最小值与最大值之间甚至有100倍以上的变化幅度。这为正确选择地下水参数带来困难，本文针对式（5）的二维地下水运动，根据推定的地下水位值，运用Guass-Newton法对大区域地下水含水层透水系数进行识别。
　　二维地下水运动的基本方程为：
　　
　　其中，h为地下水位，T为透水量系数，q为地下水涵养量或扬水量。
　　方程式（5）的透水量系数的识别，有直接法和间接法。直接法是根据已知的地下水位值，直接从式（5）中反求透水量系数T。直接法的问题是，由于地下水位中含有误差，细小的误差将导致水位的偏微分很大的误差，计算稳定性差，而且往往引起不适定问题，大大降低了参数的计算精度。本文运用间接法识别水文地质参数。即给定含水层透水量T的初始值进行迭代计算，根据3中推定的地下水位值和反复计算得到的计算水位值的残差平方和最小作为目标函数，计算最优的透水量值T。目标函数为：
　　
　　其中，hobs为ROKMT法的推定地下水位值。Tk+1为第K次迭代的透水量系数。通过下式计算，
　　
　　根据含水层的物理特性，透水量系数需满足以下限制条件：
　　
　　当透水量系数识别完成后，运用不规则三角形差分（TFDM）法进行大区域地下水计算。综合3和4的计算步骤，大区域地下水计算方法的如图2。
　　
　　5. 运用实例-Sarobetsu湿地地下水模拟
　　本文将所提出的方法运用到日本北海道Sarobetsu湿地地下水模拟中。
　　5.1 研究对象的概况
　　Sarobetsu湿地位于日本最北端（图3），为日本最大的以水台癣为主的高层湿地，面积635km2。近年来，由于农田开垦等土地开发利用，地下水位降低，湿地出现退化，海水上朔，生态环境恶化。最主要的表现是：
　　5.1.1 湿地指示性植物减少。对地下水位敏感的湿地植物水台癣面积逐渐减少，低竹类植物风长，面积逐年扩大，植物耗水量增加，植物叶面蒸散发增强。
　　5.1.2 湖泊面积减小。由于湿地干旱化导致1993年湿地内三大湖泊兜湖、Paken湖和Peken湖的面积比1975年减少10%。三湖泊周边的水台癣急剧减少，生态环境恶化。
　　5.1.3 地面沉降。湿地中心区地面最大沉降达到50厘米。
　　5.1.4 海水上朔。由于地面沉降造成和Sarobetsu河的枯水期流量减少，海水沿Sarobetsu河口上朔流入Paken湖和Peken湖，破坏了湖周边生态。
　　
　　5.1.5 Sarobetsu河下游洪水泛滥。由于Sarobetsu湿地的地面沉降和海水的顶托，造成下游洪水泛滥。
　　因此，为恢复湿地和减少地面沉降，首要课题是研究该区域地下水位恢复的措施。前提是研究该区域地下水涵养机理、地下水的运动规律。
　　5.2 地形和地质
　　Sarobetsu湿地的地形为从北到南沿海岸线逐渐降低。Sarobetsu河和海岸之间为5-20米的沙丘地带，湿地地面高程在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）。
　　地下水位的推定结果如图7。
　　10个验证点的平均误差（ME），平均绝对误差（MAE）和平方根二乘误差（RMSE）分别为0.005m，0.437m和0.598m。从三种误差看，推定的地下水位空间分布具有较好的精度，可以用来识别透水量系数。
　　5.3.3 透水系数的识别结果
　　（1）边界条件的设定
　　海岸线一侧和研究区内部的河流、湖泊作为定水头边界。其他边界位于丘陵的边沿作为流量边界。
　　（2）涵养率（入渗补给系数）和透水量系数的界限值的设定
　　
　　
　　运用分布式Tank水文模型，对1997年降雨径流的模拟，以及水平衡法水量平衡计算，降水对地下水的涵养率取0.23。根据研究区域的水文地质条件，确定透水量系数的界限值为0.0043m2/d≤T≥6800m2/d。
　　采用Guass-Newton法和有限元方法，使用非融雪期1997年6月地下水位资料，进行大区域地下水含水层透水系数进行识别。识别的透水量除以含水层厚度得到含水层透水系数，结果如图8。
　　5.3.4 98年-2000年Sarobetsu湿地地下水流动的非恒定模拟
　　为恢复湿地地下水位和生态，对1998年-2000年Sarobetsu湿地地下水流动进行了非恒定模拟。
　　（1）涵养率的设定
　　由于湿地在每年的11月底次年的5月为降雪期，融雪对地下水有很大的补给作用，本文山本等人[12]的数字积雪-融雪模型进行了模拟。模型中考虑地形标高、气温、日照量，取2℃以下时的降水为降雪，同时根据风速修正降雪量。1997年11月—2000年6月Sarobetsu湿地积雪和融雪模拟结果如图9。
　　（2）积雪和融雪的计算
　　涵养率在非积雪期取0.23，在6月至9月期间，由于植物的蒸散发强烈，涵养率取0.21。在融雪期，由于地下水主要靠融雪补给，因此，1-3月融雪量的90%作为地下水涵养，4月-5月融雪量的60%作为地下水涵养。储留系数S取0.2。
　　（3）模拟结果
　　运用TFDM法对1998年1月—2000年12月Sarobetsu湿地地下水非恒定运动进行模拟，2000年的模拟结果如图10。
　　（4）模拟结果的验证
　　对10个验证点的地下水位的非恒定计算结果进行了评价，10个点的ME为-0.063m，0.691m和0.731m。其中，第7点和第10点的误差较大，分别为1.67m和2.07m，原因是7点和第10点为丘陵山区，岩石较多，计算时涵养率设定过高。
　　
　　通过地下水运动的模拟，可以对湿地过去10-20年的地下水空间分布进行再现，结合不同年份湿地植物遥感信息和地下水空间为湿地生态恢复找到最佳地下水调控面和湿地植物退化的临界点，为湿地恢复寻求根据。其次，通过地下水运动空间分布的模拟，可以为恢复湿地的水工程的方式、规模提供依据。
　　6. 结论
　　本文在融合地质统计、逆问题理论和地下水运动理论的基础上，提出了建立大区域地下水运动模拟的理论和方法，并运用到日本Sarobetsu湿地地下水模拟中，为恢复湿地提供了科学依据。此外，本文提出的方法也可以运用到其他大区域地下水的模拟和湿地保护中。
　　参考文献
　　[1] Neuman, S.P. and E.A. Jacobson（1984）：Analysis of nonintrinsic spatial variability by residual kriging with application to regional roundwater levels, Math. Geol, 16(5), pp.499-521.
　　[2] Zhang, H.R. (1992): The abnormal problem for 2D groundwater simulation (in Chinese), J. Hydro Geo. & Eng. Geo. No.1.
　　[3] Wood, W. L. (1996): A note on how to avoid spurious oscillation in the finite element solution of the unsaturated flow equation, J. Hydrol. 176, pp. 205-218.
　　[4] Zhang, X.W., K. Takeuchi and H. Ishidaira (2001): Study on the Spurious Oscillation and Stability in Quasi Three-Dimensional Groundwater Simulation Using FEM and Comparison with SFEM and TFDM. J. Japan Soci. Hydro. Water Resour., 14(7), pp. 351-363。
　　[5] ASCE American Society of Civil Engineering Task Committee on geostatistics techniques in geohydrolgy（1990）：Review of geostatistics in geohydrology 1:Basic concepts; 2: Applications, ASCE J. Hydraul. Eng., 116(5), pp.612-658.
　　[6]Aboufirassi, M., and M.A. Marino（1983）：Kriging of water levels in the Souss aquifer, Morocco, Math. Geol., 15, pp.537-551.
　　[7] Neuman, S.P. and E.A. Jacobson（1984）：Analysis of nonintrinsic spatial variability by residual kriging with application to regional groundwater levels, Math. Geol, 16(5), pp.499-521.
　　[8] Sun, Ne-Zheng（1999）：Inverse problems in groundwater modeling, Kluwer Academic Publishers, The Netherlands, pp.161-210.
　　[9]Thomas D.（1994）：Hydroheomorphology-An introduction-,地形（日语）、第15卷A, pp.1-4．
　　[10] 张　祥伟?中津川　?竹内邦良?石平　博?山本直树?羽山早织（2002）：サロベツ湿原の非定常地下水流の解析に关する研究、北海道开発土木研究所月报、No.592、pp.13-31.
　　[11] 张祥伟、山本直树、竹内邦良、石平博、中津川诚、羽山早织（2003）情报不足条件下での广域地下水の非定常流动解析手法に关する研究，日本水文?水资源学会志、第16卷第4号、pp.349-367。．
　　[12] 山本直树?石原友海?石平　博?竹内邦良（2001）：卫星情报及び降雪?融雪モデルを用いた广域积雪水量の推定とその问题点、日本水文?水资源学会200１年研究発表会要旨集、pp.206-207．
　　作者简介：张祥伟，男，博士，毕业于日本山梨大学工学部土木环境工学科，现任中国水利水电科学研究院综合事业部副主任。Email：[email protected]
　　作者单位：中国水利水电科学研究院