趙春虎,田 干,郭國強(qiáng) ,何 淵

(煤炭科學(xué)研究總院西安研究院,陜西西安 710054)

地下水?dāng)?shù)值模擬是地下水研究中一項(xiàng)重要的應(yīng)用技術(shù),因其具有輸入數(shù)據(jù)的兼容性好,計(jì)算過程的操作性強(qiáng),輸出結(jié)果的可視化程度高等特點(diǎn),已廣泛應(yīng)用于水文地質(zhì),工程地質(zhì),環(huán)境地質(zhì)等領(lǐng)域。地下水?dāng)?shù)值模擬主要是基于有限差分法與有限元法等數(shù)值技術(shù),模型的種類有水流、水質(zhì)的預(yù)報(bào)模型、地下水資源的管理模型;模型有三維流模型、二維流模型、確定性模型和隨機(jī)性模型等。解決的問題包括地下水資源評價(jià)、地下水水位、地下水污染預(yù)測及地下水資源管理等。地下水?dāng)?shù)值計(jì)算是當(dāng)前解決地下水資源預(yù)測與評價(jià)及地下水水質(zhì)問題的強(qiáng)有力手段[1]。數(shù)值模擬作為一種定量研究地下水水量、水質(zhì)的技術(shù)手段,模型的仿真程度是影響地下水研究的主要方面,所以分析和研究地下水建模過程中誤差來源是提高模型仿真程度基礎(chǔ)工作。無論是有限差分與有限元其離散化特點(diǎn)和建模、輸入到輸出過程中不可避免存在著誤差,本文主要以當(dāng)今世界范圍內(nèi)地下水研究中應(yīng)用最廣、較為成熟的地下水模擬軟件(MODFLOW,FEFLOW等)為例,分析數(shù)值模型的誤差來源及其一般提高模型仿真程度的數(shù)值處理技術(shù)。

1 水文地質(zhì)概念模型

把具體物理模型的邊界性質(zhì)、內(nèi)部結(jié)構(gòu)、滲透性質(zhì)、水力特征和補(bǔ)給排等條件抽象化為便于進(jìn)行數(shù)學(xué)描述的概念模型,是數(shù)值模擬基礎(chǔ)工作?？陀^資料的占有情況和地質(zhì)與水文地質(zhì)條件的主觀認(rèn)識程度決定了地下水系統(tǒng)的概念模型的合理化程度。因而必須充分收集模擬區(qū)地質(zhì)與水文地質(zhì)資料,采取必要水文地質(zhì)補(bǔ)勘手段(鉆探、物探、化探、試驗(yàn)等)進(jìn)一步認(rèn)識水文地質(zhì)條件。由于含水系統(tǒng)介質(zhì)各向異性等水文地質(zhì)條件的時(shí)空差異,不可避免地進(jìn)行條件的概化,例如模型分層、邊界條件概化等。

2 數(shù)學(xué)模型

建立地下水流動(dòng)的數(shù)學(xué)模型是數(shù)值仿真技術(shù)核心內(nèi)容,主要有以水流連續(xù)性原理為基礎(chǔ)的偏微分方程與以達(dá)西定律和水均衡原理為基礎(chǔ)的水均衡方程分別構(gòu)成其數(shù)學(xué)模型,均可以得到以微分為基礎(chǔ)的微分方程(組)(式 1,以三維為例),其兩種方法結(jié)果是相同的,但后者更直觀,物理意義更明確,易于被人接受[2]。根據(jù)地下水流特征,一個(gè)由多個(gè)含水層組成的含水系統(tǒng)可相應(yīng)概化為二維模型、準(zhǔn)三維模型和三維模型[3],無疑三維數(shù)學(xué)模型更能反映地下水在三維空間中的流動(dòng)與遷移特征(壓力與物質(zhì)的遷移)。所以針對同一問題選擇不同(維數(shù))的數(shù)學(xué)模型,其誤差的原因與程度是顯而易見,實(shí)際工作中針對問題特征和不同的精度要求及現(xiàn)場資料等情況來選擇合適的數(shù)學(xué)模型是地下水?dāng)?shù)值仿真技術(shù)基礎(chǔ)。

式中：x,y,z為笛卡爾坐標(biāo)軸;t為時(shí)間;H為已知水頭;Kxx,Kyy,Kzz,為坐標(biāo)軸方向的主滲透系數(shù);μs為比彈性給水度;W為單位體積井流量,抽水時(shí)取負(fù)號;

3 數(shù)值模型的構(gòu)建

建立描述地下水流動(dòng)的數(shù)學(xué)模型后有限差分格式采用的是差商代替微商做法,源于結(jié)構(gòu)力學(xué)的有限元技術(shù)同樣通過區(qū)域剖分和插值方法將描述地下水流的定解問題化為代數(shù)方程組求解,一般采用變分原理來形成系數(shù)矩陣迭代求解,形成數(shù)值模型。所以原來數(shù)學(xué)模型的解析解(精確解),轉(zhuǎn)化為數(shù)值模型的近似解,存在著截?cái)嗾`差。這也是數(shù)值技術(shù)以離散化為基礎(chǔ)形成誤差的客觀原因。

3.1 離散化網(wǎng)格的生成問題

對研究區(qū)域在空間上的離散,平面上主要有矩形(MODFLOW、GMS等)、三角形(FEFLOW、PGMS等)、任意多邊形(陳崇希等)主要離散方式,其各具特點(diǎn)：矩形對重點(diǎn)區(qū)域如河流、斷層、抽(注)井刻畫較難,加密網(wǎng)格同時(shí)發(fā)生在行列之上,這樣相對浪費(fèi)結(jié)點(diǎn)數(shù)目,但計(jì)算格式與速度較好;三角形與任意多邊形可對重點(diǎn)區(qū)域進(jìn)行局部加密,刻畫程度較高。層間一般按含水介質(zhì)結(jié)構(gòu)概化為多層,反映在網(wǎng)格的垂直剖分上,其概化的合理程度,層間的是否反映層面空間起伏狀態(tài)均影響著模型的仿真程度。

3.2 模型的輸入中的誤差

模型輸入中的誤差來源主要是指基于數(shù)學(xué)模型的定解條件(包括邊界條件,初始條件),及模型中各項(xiàng)水文地質(zhì)參數(shù)的輸入,井孔抽注、降雨、蒸發(fā)等源匯項(xiàng)數(shù)值處理上帶來的誤差。

3.2.1 定解條件的處理

邊界條件是指模擬區(qū)邊界上的水頭分布和變化情況或者邊界上流入(出)含水層的水量分布和變化情況已知的條件。其數(shù)學(xué)形式為：

式中：H為已知水頭;Kxx為坐標(biāo)軸方向的主滲透系數(shù);μd為重力給水度;Γ1為第一類邊界;Γ2-1為潛水面邊界;Γ2-2為零流量邊界;ε′為降雨入滲補(bǔ)給量。

邊界條件主要有兩類：第一類邊界條件(已知水頭分布的邊界)(式 2,3),該類邊界條件主要技術(shù)問題在于概化的合理性,如定水頭、一般水頭邊界在時(shí)間上定量分布的輸入誤差,在空間上河流垂向切割層位、水平切割含水層的河段等輸入誤差。第二類邊界條件(已知邊界上單位寬度流量變化規(guī)律)(式 4,5),此類邊界條件在區(qū)域地下水模擬中盡量應(yīng)用其零通量邊界(隔水?dāng)鄬?、地下水分水嶺界面等自然邊界),這樣就可以極大的避免人為觀測、概化的影響,但是由于地理、地貌、地質(zhì)等條件的限制(如自然邊界距離較遠(yuǎn)),往往必須選取非零量邊界來刻畫,這樣在邊界的流量輸入的定量程度對水均衡計(jì)算有較大影響。另外在具體的研究問題中邊界條件往往存在時(shí)空不穩(wěn)定性,如河流邊界水頭季節(jié)性變化、地下水分水嶺空間推移等問題。邊界條件選擇與刻畫是影響地下水模擬仿真程度的最重要因素之一。

初始條件是指 t=0時(shí)滲流區(qū)D內(nèi)各點(diǎn)(x,y,z)處的水頭分布情況的條件。其數(shù)學(xué)形式為：

其中 H0為t=0時(shí)刻滲流區(qū) D內(nèi)各結(jié)點(diǎn)處水頭值;

初始水頭一般選取時(shí)段模擬區(qū)內(nèi)水位觀測井的統(tǒng)測水位再經(jīng)過插值(如克里格)得到全區(qū)等水頭值線作為模擬初期(t=0)初始水頭條件。在對大區(qū)域非穩(wěn)定流仿真過程中初始水頭的可靠程度極大的影響了地下水資源量評價(jià)精度,穩(wěn)定流仿真過程中初始水頭的選取影響迭代運(yùn)算速度與質(zhì)量。另外在水位統(tǒng)測過程中混合觀測井問題、動(dòng)水位問題、目的含水層水頭問題是影響水位統(tǒng)測質(zhì)量主要因素。陳崇希教授在模型應(yīng)用中采用“參數(shù)迭代法”來確定水頭未知的含水層的初始水頭提高了其合理性。

3.2.2 水文地質(zhì)參數(shù)、源匯項(xiàng)的輸入

通過野外抽水試驗(yàn)、實(shí)驗(yàn)室測試所得的水文地質(zhì)參數(shù)或經(jīng)驗(yàn)結(jié)構(gòu)參數(shù),這一系列參數(shù)是正演模型模擬計(jì)算及預(yù)測的基礎(chǔ),極大的影響管理模型、預(yù)測模型的仿真程度。在模型識別校正階段一般通過實(shí)測和計(jì)算曲線的擬合程度來調(diào)整水文地質(zhì)結(jié)構(gòu)參數(shù)(通常限定一定的取值范圍),進(jìn)一步提高模型的仿真程度。但水頭變化是地下水壓力傳遞(相對于溶質(zhì)運(yùn)移的物質(zhì)傳遞)的結(jié)果,所以其擬合標(biāo)準(zhǔn)存在一定的片面性。在模型校正過程中應(yīng)用示蹤劑示蹤(如 MODPATH),水質(zhì)模型的濃度模擬(MT3D)與實(shí)測校正等手段來綜合識別各項(xiàng)參數(shù)的合理性是更加可靠的。

源匯項(xiàng)(如降雨,蒸發(fā),開采等)其主要問題表現(xiàn)在輸入數(shù)據(jù)觀測統(tǒng)計(jì)質(zhì)量與處理技術(shù)：開采井含水層的開采層位確定,混合開采井的各含水層的流量分配,MODFLOW中主要根據(jù)導(dǎo)水系數(shù)確定各層的流量分配,黎明博士對此處理技術(shù)提出了質(zhì)疑[4],實(shí)際過程中要復(fù)雜的多。大區(qū)域開采總量的統(tǒng)計(jì)與分配,一般采用面狀分配,而不是以實(shí)際抽水井為源匯來處理的。降雨與蒸發(fā)等同樣存在處理的合理性問題,MODFLOW中降雨只發(fā)生在第一層且為瞬時(shí)補(bǔ)給地下水,一般實(shí)際過程中地下水對降雨的響應(yīng)表現(xiàn)出較明顯的滯后效應(yīng)。蒸發(fā)的處理是通過設(shè)置一經(jīng)驗(yàn)臨界深度,以此設(shè)定埋深大于此深度其蒸發(fā)量為零,這樣的數(shù)值處理方式明顯存在偏廢。陳崇希教授提出在處理混合井流和混合觀測孔時(shí),采用了“滲流 -管流耦合模型”;在處理河流入滲和降雨入滲方面采用了“入滲滯后補(bǔ)給法”;處理潛水蒸發(fā)方面,采用了非線性關(guān)系來刻畫潛水蒸發(fā)強(qiáng)度和潛水位埋深的關(guān)系,較高程度的提高了模型的仿真程度。

3.2.3 其它輸入

在我國北方地區(qū)河流是補(bǔ)給地下水的主要來源之一,類似于作為第一類自然邊界條件的河流其時(shí)空變化規(guī)律(如水位在時(shí)間上的變化,在空間上河流垂向切割層位、水平切割含水層的河段位置),在模擬區(qū)內(nèi)對河流的合理處理是相當(dāng)重要的。泉作為地下水在地表排泄點(diǎn)存在一定的技術(shù)處理難度。另外,小溪、水平擋墻、排水溝等水文地質(zhì)要素同樣存在概化的合理性問題。

4 模型的運(yùn)行與輸出

模型運(yùn)行主要體現(xiàn)在對方程組求解方法的選擇上,以有限差分為基礎(chǔ)的 MODFLOW其求解方法主要有：強(qiáng)隱式法(Strong Implic it Procedure,SIP),逐次超松弛迭代法(Successive Over Relaxation,SOR),預(yù)調(diào)共軛梯度法迭代法(Preconditioned Conjugate Gradient,PCG)[5],后者較前兩者有占內(nèi)存小,運(yùn)行速度快等優(yōu)點(diǎn),故多采用此法,但前者更適合于解決非線性問題。另外迭代次數(shù)的限制、收斂標(biāo)準(zhǔn)的設(shè)定等運(yùn)行參數(shù),對模型的計(jì)算結(jié)果必然存在影響。

5 結(jié)論

地下水?dāng)?shù)值仿真一般過程：以實(shí)地的物理模型為基礎(chǔ)經(jīng)分析概化形成合理的概念模型,選擇合適的數(shù)學(xué)模型來描述地下水流場特征,在空間和時(shí)間進(jìn)行離散化(穩(wěn)定流時(shí)間上不需要離散)形成以結(jié)點(diǎn)為基礎(chǔ)的數(shù)值模型,最后進(jìn)行模型的運(yùn)算與輸出。

從地下水系統(tǒng)的物理模型概化成合理的概念模型主要依靠于實(shí)地資料詳實(shí)程度與建模者的認(rèn)識程度與專業(yè)基礎(chǔ);選擇合適的數(shù)學(xué)模型取決于研究精度與軟件自身的限制;形成數(shù)值模型影響因素則較多,主要體現(xiàn)在模型在時(shí)空上的離散化尺度,以及各種條件輸入狀態(tài)(定解條件、參數(shù)及源匯項(xiàng)等);模型運(yùn)行可根據(jù)問題特點(diǎn)選擇不同求解方法和設(shè)定合適的運(yùn)行參數(shù)。陳崇希教授提出“防止模擬失真,提高仿真性是數(shù)值模擬的核心”[6],隨著計(jì)算機(jī)硬件水平與地下水?dāng)?shù)值仿真技術(shù)的發(fā)展,地下水模擬的仿真程度將進(jìn)一步提高。

[1]薛禹群,吳吉春.地下水?dāng)?shù)值模擬在我國——回顧與展望[J].水文地質(zhì)工程地質(zhì),1997,24(4)：21-24.

[2]薛禹群.地下水動(dòng)力學(xué)[M].北京：地質(zhì)出版社,1997.

[3]葉淑君,戴水漢 地下水流二維、準(zhǔn)三維及三維模型模擬結(jié)果比較[J]水文地質(zhì)工程地質(zhì),2003(5)：23-27.

[4]黎明,劉文波,陳崇希.MODFLOW能模擬地下水混合井流嗎?[J].水文地質(zhì)工程地質(zhì),2003,(5)：116-117.

[5]馬馳.SIP和 PCG2兩種迭代法在地下水?dāng)?shù)值計(jì)算中的應(yīng)用對比[J],西安科技學(xué)院學(xué)報(bào) 2002,22(1)：59-62.

[6]陳崇希.“防止模擬失真,提高仿真性”是數(shù)值模擬的核心[J].水文地質(zhì)工程地質(zhì),2003,(2)：1-5.

[7]Carma San Juan Kenneth E.Kolm Conceptualization,characterization and numericalmodeling of the Jackson Hole alluvialaquifer using ARC/INFO and MODFLOW Engineering Geology 42(1996)119-137.