大区域地下水模拟的理论和方法(全文).doc

大区域地下水模拟的理论和方法 张祥伟 （中国水利水电科学研究院 北京 100044） 【摘 要 （1） 其中，Lmax为准三维地下水有限元计算的三角形网格中的最小三角形最小边长；Δt为时间步长； 、分别为层间难透水层的厚度和透水系数；T1、T2为非承和承压含水层的透水量系数；μ、S分别为非承和承压含水层的储留系数。 二维的情况下，空间步长的条件为： （2） 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]，即地下水位由以下两部分构成： （3） 其中，为实测地下水位的趋势面残差， 为DEM修正后的地下水趋势面，即 （4） 其中，a0-a9为趋势方程系数，x、y为坐标，z为DEM的地形标高。 当残差满足正态性条件时，用Ordinary Kriging法进行推定，再加上式（4）即可计算得到大区域地下水的推定值。 大区域地下水位推定的方法如图1。 4. 区域地下水透水系数的识别 在地下水流动模拟中，一般根据含水层的水文地质条件和土壤特性，选定透水系数的取值。在大区域地下水计算中，含水层物理特性的信息常常十分有限，即使知道含水层的土壤类型，但由于同类土壤的透水系数的变化范围很大，最小值与最大值之间甚至有100倍以上的变化幅度。这为正确选择地下水参数带来困难，本文针对式（5）的二维地下水运动，根据推定的地下水位值，运用Guass-Newton法对大区域地下水含水层透水系数进行识别。 二维地下水运动的基本方程为： （5） 其中，h为地下水位，T为透水量系数，q为地下水涵养量或扬水量。 方程式（5）的透水量系数的识别，有直接法和间接法。直接法是根据已知的地下水位值，直接从式（5）中反求透水量系数T。直接法的问题是，由于地下水位中含有误差，细小的误差将导致水位的偏微分很大的误差，计算稳定性差，而且往往引起不适定问题，大大降低了参数的计算精度。本文运用间接法识别水文地质参数。即给定含水层透水量T的初始值进行迭代计算，根据3中推定的地下水位值和反复计算得到的计算水位值的残差平方和最小作为目标函数，计算最优的透水量值T。目标函数为： （6） 其中，hobs为ROKMT法的推定地下水位值。Tk+1为第K次迭代的透水量系数。通过下式计算， （7） 根据含水层的物理特性，透水量系数需满足以下限制条件： （8） 其中，Tjminhe和Tjmax为透水量系数的下限和上限值。 当透水量系数识别完成后，运用不规则三角形差分（TFDM）法进行大区域地下水计算。综合3和4的计算步骤，大区域地下水计算方法的如图2。 5. 运用实例-Sarobetsu湿地地下水模拟 本文将所提出的方法运用到日本北海道Sarobetsu湿地地下水模拟中。 5.1 研究对象的概况 Sarobetsu湿地位于日本最北端（图3），为日本最大的以水台癣为主的高层湿地，面积635km2。近年来，由于农田开垦等土地开发