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

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

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

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

把具體物理模型的邊界性質(zhì)、內(nèi)部結(jié)構(gòu)、滲透性質(zhì)、水力特征和補(bǔ)給排等條件抽象化為便于進(jìn)行數(shù)學(xué)描述的概念模型,是數(shù)值模擬基礎(chǔ)工作?？陀^資料的占有情況和地質(zhì)與水文地質(zhì)條件的主觀認(rèn)識(shí)程度決定了地下水系統(tǒng)的概念模型的合理化程度。因而必須充分收集模擬區(qū)地質(zhì)與水文地質(zhì)資料,采取必要水文地質(zhì)補(bǔ)勘手段(鉆探、物探、化探、試驗(yàn)等)進(jìn)一步認(rèn)識(shí)水文地質(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],無(wú)疑三維數(shù)學(xué)模型更能反映地下水在三維空間中的流動(dòng)與遷移特征(壓力與物質(zhì)的遷移)。所以針對(duì)同一問(wèn)題選擇不同(維數(shù))的數(shù)學(xué)模型,其誤差的原因與程度是顯而易見(jiàn),實(shí)際工作中針對(duì)問(wèn)題特征和不同的精度要求及現(xiàn)場(chǎng)資料等情況來(lái)選擇合適的數(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ù)號(hào);

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

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

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

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

3.2 模型的輸入中的誤差

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

3.2.1 定解條件的處理

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

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

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

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

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

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

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

3.2.3 其它輸入

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

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

模型運(yùn)行主要體現(xiàn)在對(duì)方程組求解方法的選擇上,以有限差分為基礎(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),故多采用此法,但前者更適合于解決非線性問(wèn)題。另外迭代次數(shù)的限制、收斂標(biāo)準(zhǔn)的設(shè)定等運(yùn)行參數(shù),對(duì)模型的計(jì)算結(jié)果必然存在影響。

5 結(jié)論

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

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

[1]薛禹群,吳吉春.地下水?dāng)?shù)值模擬在我國(guó)——回顧與展望[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)用對(duì)比[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.