在地下水运动数值模拟的过程中,模拟预报结果的正确与否与边界条件处理得是否恰当密切相关〔1,2〕.尤其是在人类活动影响强度日益增大的今天,在处理边界条件时,常常会面临一些新的更为复杂的问题.原因在于边界处的水流状况往往不仅受到自然因素的控制,而且还深受人类活动(如人工开采和人工补给)的影响〔3,4〕,同时还可能受到邻区水流条件变化的扰动,而对于人为边界更是如此〔5〕.所以必须对边界条件给予应有的重视,深入探讨其多重的内涵并研究出切实可行的处理方法.

　　1边界条件涵义探讨

　　在地下水运动数值模拟的过程中,一般都是在概念模型的基础上,建立描述地下水流的数学模型,然后再采用某种数值方法,对模型离散并求解.对于分布参数的地下水流数学模型而言,模型主要由两部分内容组成：①描述地下水运动规律的偏微分方程；②反映地下水模拟区域具体特征的边界条件和初始条件(若为稳定运动则没有初始条件)〔6,7〕.

　　这里的边界条件具有两重意义：一是它与初始条件一起构成地下水流数学模型的定解条件,用来说明具体目标系统的边界所具有的特定状态,从而使模型的求解能够得到切合实际状况的特解.二是它描述了目标系统与其周围环境之间的相互作用关系,即它们之间物质、能量和信息的交换.也就是说,周围环境的变化会对目标系统产生影响,而目标系统的变化,也会对周围环境产生作用.

　　边界条件的第二重意义的重要性伴随着水资源匮乏的日益严重而逐渐凸显出来.一方面,在过去,由于开采规模往往比较小,人工开采对于边界水流状态的影响还比较小,可以忽略不计.因此,在边界处可近似将水流视为仅受自然因素影响,这样的边界条件处理起来比较简便.然而,随着开采规模的扩大,人工开采对边界处的水流已造成强烈影响,使边界处的水流状态同时受到自然、人为双重因素的制约.另一方面,在过去,往往周围邻区尚不存在其它的开采水源地.而现今面临的情况是：在相邻区域,在目标系统的上游或下游存在其它开采水源地,相邻水源地常常各自划定边界并分别建立模型,它们在边界处的耦合效应要通过边界条件给予描述〔8,9〕.

　　在地下水运动数值模拟过程中,边界条件的处理可划分为两个阶段.第一个阶段即在模型识别检验阶段,边界条件利用以往的观测资料即可确定,处理起来较为容易.第二阶段即在模型预报阶段,边界条件处理起来较为复杂.这里存在首先要对边界条件做出预报的问题,然后才能进行模型预报.因为在模型预报时,边界条件是作为已知条件出现的.边界条件的预报既要考虑本区和邻区未来自然条件的变化,同时也要考虑本区和邻区未来人为影响(如人工开采和人工补给)的变化,即边界条件的预报要综合考虑各种作用(包括自然因素、人为因素及相邻系统的扰动等)的复合效应.否则,边界条件就不能刻画边界处水流状态的真实变化,并将最终导致整个模型计算结果的不正确.

　　2处理方法实例

　　下面给出两种处理边界条件的方法和实例,它们分别都刻画出了自然因素和人为作用对于边界水流的复合影响.

　　2.1北京市平谷电厂水源地数值模拟一类边界水位的预报在该水源地的集中人工开采没有上马投产之前,边界处的水位只受降水、地质条件等自然因素控制,通过建立回归方程

　　h末=α0+α1P+α2h初

　　来表达某时段的末刻水位h末与该时段内的降水量P和初始水位h初之间的数量关系.其中α0为常数,α1、α2分别为回归系数.这里h末只受自然因素控制.在该水源地的集中人工开采投产之后,集中人工开采的影响,这里主要是指由于大量人工抽水而引起的水位降落漏斗,要扩展到离抽水井较近的一类边界所在的位置,其最近距离为r=10km,按设计开采量Q=80000m3／d,可根据?Theis?公式近似算出离开采井最近的一类边界处的水位降深变化状况如表1所示.

　　由此可以看出,一类边界处的水位不仅受自然因素控制,还要受到集中人工开采的影响,是二者叠加复合作用的结果.在预报阶段,在模型预报之前首先对边界条件做出预报.作者采用的方法是根据预报阶段的设计开采量,运用?Theis?公式近似算出集中人工开采在一类边界上引起的水位降深S,然后运用下式计算某时段末刻水位h?末(自然影响)与水位降深S(人为影响)的代数和(二者叠加),从而确定H末：

　　H末=h末-S

　　H末即为在h末的基础上叠加了人工影响后该时段的末刻水位,该水位综合反映了自然因素和集中人工开采对一类边界水位的共同影响.

　　2.2某水源地数值模拟二类边界单宽流量的预报该水源地的含水层为我国北方中、下奥陶统的碳酸岩岩层.南部山区灰岩裸露,地势较高,接受降水补给,为岩溶潜水.北部丘陵区地势较低,灰岩含水层为石炭系和二叠系砂页岩岩层所覆盖,向北埋深逐渐增大,为岩溶承压水.地下水总体流向为由南向北.在划定模拟区域时,在北部边界上,由于承压岩溶含水层深伏地下,根据实际情况将其确定为二类弱透水的排泄边界.同样在该水源地没有上马抽水之前,二类边界处的单宽流量q只受降水、地质条件等自然因素控制.而在水源地上马抽水之后,由于人工抽水而产生的水位降落漏斗会波及到离抽水井较近的二类边界所在的位置,其最近的距离为r=1375m,按设计开采量Q=10000m3/d,根据?Theis?公式近似算出离开采井最近的二类边界处的水位降深变化状况如表2所示.

　　表1距离开采井最近一类边界处的降深变化

　　抽水时间／d60100365730

　　水位降深／m0.0050.0180.1180.203

　　表2距离开采井最近的二类边界处的降深变化

　　抽水时间／d3060100365

　　水位降深／m0.4450.5000.5410.645

　　可见,此处边界上水位的实际变化,以及由此诱导出来的水力梯度I和单宽流量q的实际变化,也是自然因素与人为作用复合影响的结果.在模型预报之前,也必须首先要对边界条件进行预报.对该水源地北部的二类边界而言,就是要对单宽流量q做出预报.考虑到在北部边界处承压含水层深伏地下,可用资料很少,作者采用迭代逼近方法来预报北部二类边界的单宽流量.

　　由达西定律,单宽流量q=T·I,其中T为导水系数,I为水力梯度.由于北部二类边界位于承压区,边界处的承压含水层厚度M和渗透系数K均为定值,故T=KM可视为常数.因此,单宽流量仅随水力梯度I的变化而变化.无论是天然降水还是集中人工开采状况的变化,它们对北部二类边界单宽流量的作用都最终表现为对边界处水力梯度的改变.所以,在预报时只要抓住水力梯度这个关键因素即可.

　　具体步骤是：对于第i个区段,先给出水力梯度的初估值,由算出相应的单宽流量,据此运转预报模型,并根据预报结果(水位值)计算出北部二类边界附近的水力梯度,记为,然后再运用根据公式算出,并再运转预报模型.重复进行上述迭代过程,直到前后相邻两次迭代中的单宽流量的差值的绝对值小于预先给定的允许误差ε为止,即要满足条件：