龔嘉臨， 朱 松， 李非里， 金贊芳

(1．浙江工業(yè)大學(xué)環(huán)境學(xué)院， 浙江 杭州 310032；2.浙江省水利水電勘測設(shè)計院， 浙江 杭州 310002)

0 引言

目前，全球幾乎所有湖泊面臨著富營養(yǎng)化及其引起的藍藻水華問題，嚴重影響著湖泊的水質(zhì)和生態(tài)功能[1-2]。在國內(nèi)，湖泊富營養(yǎng)化問題也早已成為一個較為突出的環(huán)境問題，如太湖、滇池、鄱陽湖等大型湖泊均已處于富營養(yǎng)或重富營養(yǎng)狀態(tài)。而且，無錫飲用水供水危機等事件的發(fā)生也表明了我國湖泊水體富營養(yǎng)化問題的嚴重。然而，盡管發(fā)達國家對湖泊的富營養(yǎng)化研究和治理已經(jīng)經(jīng)歷了半個世紀，但污染湖泊的水質(zhì)和生態(tài)恢復(fù)至今仍未能完全實現(xiàn)。這是因為湖泊的富營養(yǎng)化過程是一個不可逆且長期產(chǎn)生影響的過程[3]。首先，造成湖泊富營養(yǎng)化的外源有很多種，包括農(nóng)業(yè)污染、工業(yè)污染、生活污水，大氣沉降、地下水入湖等。其次，湖泊底泥和腐爛的水生植物等內(nèi)源物質(zhì)也會將富集的營養(yǎng)物質(zhì)釋放，從而影響湖泊。而在地表水源受到普遍污染的現(xiàn)在，人們進一步開發(fā)并利用地下水，導(dǎo)致淺層地下水逐漸成為氮元素的貯存庫。其中，硝酸鹽污染更是已成為地下水污染的重要問題[4-5]。由于地下水與湖泊緊密聯(lián)系，存在著水力交換，所以無論是哪種系統(tǒng)發(fā)生變化，都會對關(guān)聯(lián)的另一系統(tǒng)造成影響。

隨著對湖泊生態(tài)系統(tǒng)的不斷研究，湖泊周邊地下水對湖泊富營養(yǎng)化問題的影響受到了國內(nèi)外研究人員的廣泛關(guān)注。自上世紀70年代，國外的研究人員就開始了對地下水與湖泊之間相互作用關(guān)系的研究。由于賦存于地面下的巖石空隙中的地下水在一般情況下不能被人們直接觀察，所以研究人員只有通過在現(xiàn)場的水文地質(zhì)勘察與地下水動態(tài)的持續(xù)監(jiān)測才能了解區(qū)域地下水系統(tǒng)的存在形式及運動規(guī)律。在研究初期，研究人員對地下水-湖泊系統(tǒng)的研究手段主要是利用滲流計和壓力計等現(xiàn)場監(jiān)測儀器來估算淺層地下水和湖泊之間的交換量。

隨著計算機技術(shù)的發(fā)展和利用，越來越多的研究人員使用可視化模型來定量化模擬區(qū)域地下水和關(guān)聯(lián)湖泊之間的相互作用[6-8]。其中，有限差分模型MODFLOW作為應(yīng)用最廣泛的模型軟件，受到了眾多研究人員的青睞。本文通過對地下水水流模型MODFLOW、溶質(zhì)運移模型MT3DMS和湖泊模擬方法的介紹，總結(jié)其功能與應(yīng)用特點。

1 有限差分模型MODFLOW簡介及特點

MODFLOW模型（Modular Three-dimensional Finitedifference Ground-water flow model）最初是由美國地質(zhì)調(diào)查局（US Geological Survey）的MCDONALD等在上世紀80年代所開發(fā)出來的一款數(shù)值模擬程序，用于模擬孔隙介質(zhì)中地下水流動。

MODFLOW基于離散化方法，將連續(xù)的區(qū)域空間和時間分為一系列離散的格點。在一定精度下，將區(qū)域空間離散為有限個規(guī)則單元的集合體，使得每個網(wǎng)格單元上的各種參數(shù)值近似為常數(shù)。在離散化之后，整個區(qū)域的計算問題就被概化為了有限個單元網(wǎng)格的計算。結(jié)合時間離散之后，利用有限差分方程構(gòu)成線性方程組。對線性方程組聯(lián)立求解之后，獲得各個離散點上的水頭數(shù)值近似解。這種方法可以較好地反映區(qū)域地下水的狀態(tài)，且具有較高仿真度。

MODFLOW的一個最顯著特點是其模塊化的結(jié)構(gòu)[9]。它在主程序的基礎(chǔ)上可以整合一系列相對獨立的子程序包，而這些子程序包是由具有類似功能的子程序所組成的。同時，研究人員可以根據(jù)實際情況選用其中合適的子程序包對地下水系統(tǒng)進行數(shù)值模擬。這種模塊化結(jié)構(gòu)讓整個程序更容易理解、修改，并可對其添加新的子程序包，如河流子程序、模擬水平流動障礙的子程序包等。模塊化的結(jié)構(gòu)使得MODFLOW的功能和應(yīng)用范圍大大拓展。模塊化結(jié)構(gòu)的MODFLOW不僅僅局限于最初的模擬地下水在孔隙介質(zhì)中的流動，在合理使用的情況下可以解決裂隙水的流動等實際問題。

目前，使用有限差分模型MODFLOW的軟件主要有GMS（楊百翰大學(xué)環(huán)境模型實驗室和美國軍工水道實驗室合作開發(fā)），Visual Modflow（加拿大Waterloo水文地質(zhì)公司開發(fā)）和Visual Modflow flex。隨著地下水數(shù)值模擬軟件的廣泛應(yīng)用，國內(nèi)外研究人員主要通過MODFLOW模擬和預(yù)測區(qū)域地下水環(huán)境的動態(tài)變化，實現(xiàn)對實際問題的評價與動態(tài)預(yù)測，從而制定合理的措施解決問題[10]。SONG等[11]通過構(gòu)建模型來確認開采計劃是否合理，從而確保巴拉素水源地的安全。馬馳等[12]根據(jù)模型結(jié)果預(yù)報西華水源地地下水位，獲得合理的水源地開采量，為西北地區(qū)地下水資源的評價開創(chuàng)了先例。尉鵬翔[13]利用Visual Modflow軟件對北京某污染場區(qū)進行建模（利用MODFLOW和MT3D），分析了污染物的遷移，為該污染場區(qū)地下水評價和污染治理提供了科學(xué)依據(jù)。ELKRAIL[14]結(jié)合GIS和MODFLOW對地下水的脆弱性進行評價，用于指導(dǎo)松花江所在區(qū)域地下水資源的開發(fā)利用和保護。

因為MODFLOW的最終目的是為了模擬實際的水文地質(zhì)情況，故在使用該模型之前需要盡可能多地收集研究區(qū)的水文地質(zhì)資料。在此基礎(chǔ)上，整個模型的構(gòu)建步驟見圖1。

圖1 MODFLOW模型構(gòu)建一般步驟

2 溶質(zhì)運移模型MT3DMS

為了研究地下水和湖泊富營養(yǎng)化的關(guān)系，在利用MODFLOW獲得水流模型的基礎(chǔ)上，研究人員還需要根據(jù)溶質(zhì)運移模型獲得相應(yīng)的水質(zhì)模型來進一步估算地下水對湖泊的污染負荷。

20世紀90年代，ZHENG開發(fā)了模塊化的三維溶質(zhì)運移模型MT3D。這款模型開發(fā)完成之后就廣泛地應(yīng)用于地下水中的溶質(zhì)運移模擬和研究。隨后，可插入多重污染物組分生化反應(yīng)模塊程序的MT3DMS由ZHENG等[15]于1998年開發(fā)。MT3DMS的功能在MT3D的基礎(chǔ)上得到了極大的擴散，其獨特之處在于囊括了有限差分法、基于粒子示蹤的歐拉-拉格朗日法以及高階有限差分TVD法3類主要的溶質(zhì)運移解算法。之所以將3種解算法綜合起來是因為沒有一種單獨的方法可以使用所有的運移條件[16]。

與MT3D一致，MT3DMS也是以塊體中心法有限差分水流模型（如MODFLOW）為背景開發(fā)的。同時，MT3DMS保留了MT3D的模塊化結(jié)構(gòu)。而且和MODFLOW一樣的是，MT3DMS同樣是模塊結(jié)構(gòu)。其模塊化的結(jié)構(gòu)可以單獨模擬運移的對流、彌散/擴散、源匯問題和化學(xué)反應(yīng)問題，有利于合理分配內(nèi)存；處理其它運移過程與化學(xué)反應(yīng)的模塊可以隨時插入程序，而不需要改動現(xiàn)有代碼。因此，MT3DMS可以模擬含水層中的對流、彌散-擴散過程，也可以模擬地下水污染物質(zhì)在運移過程工程的生物反應(yīng)和化學(xué)反應(yīng)。趙貝等[17]通過MODFLOW和MT3DMS構(gòu)建了河北省東南部地區(qū)的一個垃圾填埋場的三維數(shù)學(xué)模型，并模擬預(yù)測了該區(qū)未來30 a內(nèi)的特征污染物Cl-的遷移情況。模擬結(jié)果表明，污染物濃度隨著離污染源的距離增大而減少，主要危害地下水流方向上的地表水體和淺層地下水；湖泊對污染物具有持續(xù)捕獲作用；開采地下水會影響污染物的遷移速率和方向。GUSYEV等[18]利用MODFLOW，MT3DMS和同位素氚對新西蘭的西陶波湖捕集區(qū)進行地下水年齡模擬，用于構(gòu)建流域地下水年齡分布。

3 湖泊模擬方法

自然環(huán)境中的地下水和湖泊存在著水力交換，從而導(dǎo)致地下水影響著湖泊的水平衡、營養(yǎng)化水平和酸緩沖能力。而湖泊周邊的地下水是作為補給源還是排泄匯，這取決于湖泊與周邊地下水的水力梯度。1979年，BORN等[19]依據(jù)湖泊流域的水文循環(huán)特征，對地下水-湖泊交換定義了3種不同的類型，見圖2。第1種類型是地下水穿流，湖泊的部分湖底接受地下水，而部分湖底發(fā)生滲漏；第2種類型是地下水進入湖泊，對湖泊進行補給；第3種是湖泊水對地下水進行補給。一般情況下，地下水-湖泊交換類型隨季節(jié)變化而發(fā)生變化。其中，第1種類型在山區(qū)流域中出現(xiàn)的頻率較多，且在具有恒定流向地表水體的平原地區(qū)也常出現(xiàn)，具體表現(xiàn)為湖泊在上游的入湖區(qū)接受地下水和下游的出湖區(qū)滲漏進入地下水；第2種類型常出現(xiàn)在沙漠地區(qū)或蒸發(fā)量遠遠大于降雨量的湖泊流域，當湖泊的周邊地勢較高時，湖泊位于流域的地勢較低處，湖泊水體的出流主要依靠蒸發(fā)的作用，入流則主要來自部分降水和地下水流；第3種類型通常在地表水位主要補給的湖泊流域出現(xiàn)，湖泊水體的出流受地勢影響，或湖底的滲透性能較好，從湖底補給的地下水能夠較快地排走，從而使含水層中的地下水位維持在比湖泊水位更低的水平。

圖2 地下水-湖泊的交換形式

研究地下水和湖泊的相互作用有利于了解整個系統(tǒng)，從而為解決實際問題提供依據(jù)。而在數(shù)值模擬中，研究人員也可以通過以下幾種方式模擬地下水和湖泊的相互作用[20-22]。

3.1 湖泊區(qū)域作為高滲透系數(shù)的含水層

當進行MODFLOW模擬時，將湖泊作為構(gòu)建的網(wǎng)格模型含水層的一部分，同時設(shè)置其所在區(qū)域的滲透系數(shù)為高滲透系數(shù)。李勇[23]和莫美仙[24]結(jié)合該方法和第2種方法，把滇池作為模型含水層中的一部分，設(shè)置湖泊所在區(qū)域單元具有高滲透系數(shù)，且將湖泊區(qū)域設(shè)置為邊界，從而來模擬滇池。前者利用ZONGBUDGET模塊計算出滇池東側(cè)地下水入湖水量為8 073.8 m3/d，硝態(tài)氮和總磷入湖量分別為32.4和0.91 t/a，北側(cè)地下水入湖水量為8 463.6 m3/d，硝態(tài)氮和總磷入湖量分別為54.72和0.832 t/a；后者計算出昆明盆地北部地下水入湖量在6.5×106～6.9×106m3/a之間，并根據(jù)地下水攜帶的氮、磷濃度平均值乘以入湖水量估算氮和磷的入湖量范圍分別為27.13～ 406.26，1.84～ 10.04 t/a。

這種方法模擬較為便利，可以計算湖泊水位，并模擬穩(wěn)態(tài)和瞬態(tài)問題。但其缺陷在于只能模擬滲漏湖泊，無法包含其它的地表水流動，且容易出現(xiàn)質(zhì)量不守恒問題。

3.2 湖泊作為邊界條件輸入模型

湖泊作為邊界條件輸入模型即把湖泊這一系統(tǒng)看成是MODFLOW模擬區(qū)域的其中一個已知邊界。在這種情況下，研究人員可以利用MODFLOW中的定水頭邊界或通用水頭邊界來表示湖泊。潘田[25]將研究區(qū)東北部連接的太湖定義為定水頭邊界，以太湖水位作為水頭條件，從而實現(xiàn)研究區(qū)地下水與太湖的水力聯(lián)系。同時，結(jié)合MT3DMS進行溶質(zhì)運移模擬獲得太湖流域地下水硝態(tài)氮補給太湖。結(jié)果表明，太湖周邊淺層地下水補給太湖水量為0.2億m3/a，占總排泄量的6%；太湖周邊淺層地下水中的硝酸鹽氮入湖負荷約為14 245 t/a，占補給太湖氮污染負荷的27.8%。郭玉川[26]將博斯騰湖作為變水頭邊界，湖泊水頭值根據(jù)實際水位分時段賦值，通過構(gòu)建模型，從而分析模擬了基于生態(tài)安全水位的不同調(diào)控方案的水均衡情況，對當?shù)氐叵滤坏恼{(diào)控方案進行優(yōu)選。

該方法模擬的湖泊水位穩(wěn)定且利于模擬滲漏和排水湖泊。但是，由于湖泊水位已經(jīng)作為固定邊界輸入，在整個模型模擬過程中無法體現(xiàn)湖泊的水位變化且一般不包括其它地表水流。同時，當其作為定水頭邊界輸入模型的時候，湖泊邊界其實相當于無限的源匯項。

3.3 河流程序包代替湖泊模塊

河流模塊作為MODFLOW邊界條件中的常用模塊，也可以用于模擬其它地表水體。這種方法與通用水頭邊界條件一樣都是以外部應(yīng)力形式導(dǎo)入模型。但與通用水頭邊界導(dǎo)入不同的是，使用河流模塊模擬湖泊等其它地表水體的時候可以定義水體范圍，同時在結(jié)果中可以顯示水域的變化情況及流場。而缺陷在于可能出現(xiàn)質(zhì)量不守恒情況。

3.4 LAK3程序包

上文已說明MODFLOW是一種模塊化的結(jié)構(gòu)，這為湖泊程序包的開發(fā)和使用提供了便利。LAK湖泊程序包（LAK1,LAK2,LAK3）從上世紀90年代開發(fā)至今，已被國內(nèi)外研究人員廣泛應(yīng)用。KROHELSKI等[27]利用 MODFLOW 和LAK3程序包來模擬威斯康星州西北部Dane鎮(zhèn)的3個淺水滲透湖泊（Fish湖，Mud湖和Crystal湖）及其周圍的淺層地下水。模擬結(jié)果表明，F(xiàn)ish湖水位在抽水條件下會降低；在沒有持續(xù)抽水的情況下，湖泊水位又會在幾年內(nèi)恢復(fù)。FEINSTEIN等[28]模擬了威斯康星州西北部Waukesha鎮(zhèn)的湖泊、井和淺層地下水的相互作用。模擬結(jié)果表明，地下水是Beaver湖最重要的入流補給部分（等于總流入的59%）；對于Pine湖和North湖來說，地下水仍然是重要的來源（分別占總流入的16%和5%），但小于降水和地表水的貢獻。當在地表61 m以下以213.85 L/min速率抽取地下水5 a，對湖泊水位的影響也并不大，模擬的Pine湖的水位和出流量僅分別減少了0.91 cm和3%。但是，加入抽水速率增加到910 L/min，模擬的Pine湖的出流量明顯減少，相對基本流量條件下降了14%。在我國，應(yīng)用LAK3模塊的情況較少。程春龍等[29]利用MODFLOW和LAK3模塊建立了天津擬建寧車沾水庫的耦合模型。模擬結(jié)果表明，水庫水位下降速率受到當?shù)氐恼舭l(fā)強度影響較大，受滲漏影響極小，年水位降幅1.16 m。

LAK3程序包作為最新的模擬湖泊程序包，具有進行穩(wěn)態(tài)和瞬態(tài)模擬的功能，且可以用于模擬地表水流和溶質(zhì)運移。同時，和河流程序包一樣，LAK3程序包的湖泊模擬可以在結(jié)果中顯示。與前面3種方法相比，LAK3輸入較為復(fù)雜，包括湖泊水位、底泥層厚度、滲透系數(shù)、降雨、蒸散、徑流等參數(shù)。其缺陷主要有2點：①當滲透系數(shù)較高時，模擬具有潛在的不穩(wěn)定性；②在沒有用戶干涉的情況下，可能出現(xiàn)錯誤結(jié)果。

4 存在的問題及發(fā)展趨勢

根據(jù)文獻檢索發(fā)現(xiàn)，近年來有關(guān)區(qū)域地下水-湖泊系統(tǒng)的三維可視化模型研究數(shù)并不少。這些研究涉及的湖泊及周邊地下水的水文地質(zhì)條件復(fù)雜多樣。利用MODFLOW等模型的構(gòu)建為區(qū)域地下水和湖泊的水資源預(yù)測規(guī)劃和管理提供了重要的依據(jù)。但是，綜合分析可以發(fā)現(xiàn)以下問題的存在。

（1）研究范圍常以人為邊界劃定。由于實際地質(zhì)條件、資料收集、工作量和計算機容量所限，研究人員在進行模擬研究的過程中往往很難找到自然邊界作為邊界條件，這就導(dǎo)致模擬區(qū)域往往比實際的地下水含水層系統(tǒng)要小。這種情況下，MODFLOW模型常采用人為邊界。因為人為邊界的引入，也導(dǎo)致了模型預(yù)測的困難。

（2）水文數(shù)據(jù)資料收集和驗證困難。和國外相比，國內(nèi)的野外實驗基地沒有形成規(guī)模，這也導(dǎo)致水文地質(zhì)資料往往沒有長期的監(jiān)測、考證資料，無法形成數(shù)據(jù)積累。這也導(dǎo)致了國內(nèi)構(gòu)建的模型難以進一步修正優(yōu)化，達不到長期預(yù)測的效果。

（3）數(shù)據(jù)處理往往以多年平均、短期觀測資料為主。同樣是因為實際地質(zhì)條件、資料收集和工作量所限，研究人員往往只能盡可能將收集到的水文地質(zhì)信息進行合理概化，以求在模擬的時候接近實際情況。因此，這也造成了根據(jù)數(shù)據(jù)資料的差異，不同研究人員的研究成果往往存在一定差異。同時，這也必然導(dǎo)致研究成果與實際的偏差，以至于很多情況僅僅限于科學(xué)研究。

5 結(jié)論

MODFLOW這類用于模擬實際環(huán)境狀態(tài)的模型離不開各種大數(shù)據(jù)的支持。隨著我國社會經(jīng)濟的不斷發(fā)展，對于城市及周邊地區(qū)地質(zhì)數(shù)據(jù)的調(diào)查和完善提出了新的要求。進入21世紀以來，中國地質(zhì)調(diào)查局和各地方政府合作開展的“城市地質(zhì)調(diào)查項目”就是為了解決城市發(fā)展過程中的地質(zhì)問題，并已在多個重要城市展開試點調(diào)查工作。這為將來的進一步地質(zhì)研究提供了便利。而隨著互聯(lián)網(wǎng)的發(fā)展，大數(shù)據(jù)資料庫的建立使區(qū)域水文地質(zhì)資料得到了不斷完善和擴充，加快了信息交流。同時，計算機技術(shù)的發(fā)展也使計算機容量不斷增大，為更大范圍的區(qū)域地下水模擬提供了可能。長遠來看，數(shù)據(jù)資料的擴充和共享交流對模型研究具有重要意義。