高放廢物深地質處置地下水流數值模擬方法研究進展

李露露，周志超，邵景力，崔亞莉，趙敬波

（1.中國地質大學(北京)水資源與環境學院，北京 100083；2.核工業北京地質研究院，北京 100029）

核能已成為人類社會中必不可少的能源之一。然而，核科技與核能技術的發展伴隨產生了一定的放射性廢物，特別是高水平放射性廢物（簡稱高放廢物），具有放射性極強、毒性極大、持續時間長等特點[1]，一旦進入人類生存環境，危害極大且難以消除，因此，如何對其安全處置仍然是世界性難題。目前國際上公認的最主要的處置方式是深地質處置，是指將固體形式的高放廢物埋藏在距地表500～1 000 m 的地質體中，使之永久與人類生存環境隔離[2]。然而，在放射性廢物最終處置過程中，有害的放射性核素遷移到人類環境中最有可能的方式就是地下水的搬運。處置場建成運轉后，填埋廢物體的核素釋放到地下水中，并隨之遷移擴散，造成地下水甚至地表水體的污染。因此，要有效地進行高放廢物的深地質處置，就要了解區域深部地下水形成演化機制，利用數值模擬技術預測未來相當長時間系統內地下水的流動特征已逐漸成為高放廢物處置庫選址和場址評價中必不可少的研究內容。

放射性核素在地質體中的遷移機理及過程極為復雜，有對流、彌散、吸附、衰變和微生物等物理、化學反應過程，目前準確模擬核素在介質遷移的全過程仍然十分困難。以上過程中，核素隨地下水流遷移是最主要、最直接的方式，也是分析和了解其他過程的基礎。由于處置場巖石空隙類型、構造以及對地質和水文地質資料掌握程度的不同，建立的水文地質概念模型和數學模型亦不相同。此外，由于深部資料相對匱乏，加之基巖含水介質高度的非均質和各向異性，地下水流數值模擬結果會出現較大的不確定性，因而不確定性分析工作必不可少。現階段已有很多可用于地下水流數值模擬的軟件，如何選擇適用的軟件也是值得關注的問題。

本文就以上問題展開綜述，以了解世界各國高放廢物深地質處置地下水流數值模擬方法研究進展及其應用，可為我國高放廢物處置場選址、安全評價和后續放射性核素遷移轉化的預測模擬提供基礎支撐。

1 深地質處置地下水流數值模擬方法

深地質處置中涉及的圍巖類型大多為低滲透性基巖，常見的有花崗巖、黏土巖、凝灰巖和巖鹽等。巖體中的含水空間主要有裂隙介質或裂隙-孔隙雙重介質。孔隙發育較為均勻，而裂隙的隙寬、展布方向、延伸長度等均具強烈的非均質性和各向異性。因此，深入了解裂隙巖體滲流理論及過程是深地質地下水流數值模擬的基礎。目前用于描述裂隙巖體滲流的模型可分為等效連續介質模型、離散裂隙網絡模型、雙重介質模型、等效-離散耦合模型[3]。

1.1 等效連續介質（EC）模型

EC 模型是將裂隙中的水流等效平均到整個巖體中，將裂隙巖體看作具有對稱滲透張量的各向異性連續體，利用經典的、成熟的連續介質理論進行分析。實際上就是求解地下水在多孔介質中流動的控制方程：

式中：H—水頭/m；

K—滲透系數張量/（m·d?1）；

μs—貯水率/m?1；

W—源匯項/d?1。

這類模型適用于大區域、長序列、裂隙發育程度較高或較均勻的地區，其優點是技術方法成熟、需要數據獲取較為容易、操作性好。但等效連續介質方法存在著明顯的缺陷：不能很好地刻畫裂隙的特殊導水作用，同時由于裂隙的非連續特征，不是所有的裂隙巖體都可以等效成連續介質，典型單元體的大小和等效水力參數較難確定。因而適用范圍有限，不適合場地和儲罐區的裂隙水流的精細模擬。Long 等[4]認為應用EC 模型首先要確定裂隙巖體的等效滲透張量，指出滲透張量必須要無條件地適用于動力場相似的水流系統。Long 等[5]進一步考慮采用幾何形態法確定滲透張量，探討了裂隙間連通程度對裂隙滲透率的大小和性質的影響。Snow[6]通過假定裂隙面是無限延伸的，推導了單個裂隙、一組平行裂隙的滲透張量公式，并提出多組裂隙巖體的滲透張量可以疊加。Oda[7]將包含大量地質不連續面的巖體視為各向異性、彈性多孔介質，建立了求解應力和流體流動耦合的控制方程組，并利用現場實例驗證了該模型的適用性。田開銘等[8]建立滲透張量的改進模型，并將統計測量與現場試驗方法求得的滲透張量進行對照和校正，提高了滲透張量計算的準確性和可靠性。

1.2 離散裂隙網絡（DFN）模型

DFN 模型是以裂隙網絡中的單個裂隙為研究對象，認為基質巖塊本身不透水，地下水只在裂隙中流動，將裂隙網絡視為非連續體。以單個裂隙內水流基本公式為基礎，利用流入和流出各裂隙交叉點流量相等的原則建立方程，通過求解方程組獲得各裂隙交叉點的水頭值。按照Wittke 線素模型[9]的基本原理，把裂隙交叉處作為節點，節點之間的裂隙稱為線單元，每個線單元流向共同節點的流量等于0（穩定流）或等于儲存量的變化量（非穩定流），因此表征裂隙巖體滲流網絡中i節點單元域內地下水均衡方程式為：

式中：qj(j=1,2,···,N′)— 表征單元域內某一時刻流進流出各銜接線元的流量/（m3·d?1）；

wj(j=1,2,···,N′)—表征單元域中每個線元上的垂向補給量/（m3·d?1）；

Qi—i節點上的源匯項/（m3·d?1）；

Hfi—節點上的水頭/m；

N′—線單元個數；

Si—裂隙以i點為中心的表征單元域內彈性貯水（釋水）系數；

ej—裂隙的隙寬/m；

lj—線元的長度/m。

DFN 模型比較真實地描述了裂隙網絡中地下水流特性，比較適合解決中小尺度（處置場地、儲罐尺度）、需要精細刻畫的地下水流問題。Wilson 等[10]提出了兩種模擬二維裂隙網絡水流的有限元技術，指出當裂隙交叉點處的水流干擾小到可以忽略不計時，可以使用線性單元程序。三角形單元程序可用于確定基質-裂隙耦合流。Wittke 等[9]提出了類似電路分析回路中的網絡線素法來模擬裂隙網絡。王明玉等[11]利用有限元方法建立三維多邊形DFN 模型，模擬裂隙巖體非均質性和各向異性等滲流特征，進而確定裂隙巖體典型單元體及滲透張量的大小。Long 等[12]提出了三維圓盤隨機DFN 模型，假定裂隙為不可滲透基質中的圓盤狀不連續面，使用混合解析-數值技術計算通過該裂隙網絡模型的穩定流量。Nordqvist 等[13]提出了可變隙寬的DFN 模型，用該模型研究了裂隙巖體中水流和溶質運移規律。但大多數情況下，盡管在處置場區、儲罐區投入大量的勘探、測量、試驗等相關工作，但很難準確獲得每條裂隙的裂隙特征（裂隙的產狀、張開度、間距、跡長等）及其連通性，同時存在工作量大、耗時多的缺點。

到目前為止，復雜的三維裂隙網絡滲流，尚無成熟的算法和模擬軟件，通常利用隨機模擬生成裂隙網絡并進行滲流模擬。該方法首先要形成DFN，但由于裂隙系統的非均質各向異性，裂隙的分布和特性存在著很大的不確定性，因而發展出統計學方法生成DFN。目前常用的統計學方法是Monte-Carlo 法，通過實測地表裂隙數據的統計規律反推地下空間巖體裂隙幾何參數近似解。DFN 的建立首先需要測量和收集裂隙特征空間分布數據，形成裂隙面幾何參數的概率統計模型，應用Monte-Carlo 隨機方法生成DFN，最終對模型構建結果進行驗證，調整模型參數確定更接近實際情形的DFN[14]。Andersson 等[15]運用三維裂隙網絡模型模擬了不同水力條件下的滲流過程。Dverstorp 等[16]建立了稀疏裂隙巖石中的DFN 模型，并模擬了該巖石中的示蹤試驗結果。何楊等[17]針對稀疏的三維主干巖體裂隙網絡，以滲流區內的水量平衡原理為依據推導了三維巖體裂隙網絡非穩定滲流的基本方程，并利用Monte-Carlo 模擬了巖體中的裂隙網絡。宗自華等[18]基于已有地表裂隙統計數據，對北山BS03 鉆孔周邊巖體進行了裂隙網絡建模，并分析了其連通性。黃帆等[19]基于Monte-Carlo 法，引入裂隙跡線長度和開度的相關性函數，提出改進的隨機裂隙網絡生成算法，驗證了算法的有效性和程序的正確性，最終借助DFN滲流模型對生成的裂隙網絡進行滲流分析。

1.3 雙重介質模型

雙重介質模型認為裂隙巖石是由大裂隙系統和巖塊孔隙或細小裂隙系統共同構成的連續介質的統一體，大裂隙導水性強，主要起到地下水流的通道作用，稱為主干裂隙介質。巖塊孔隙或細小裂隙分布均勻儲水空間大，主要起到儲存和釋放地下水的作用，稱為孔隙介質或分枝裂隙介質。根據實際問題的特點以及模擬時的要求，又分為：（1）雙重孔隙度模型。該模型認為很多情況下，巖塊的滲透率遠小于裂隙的滲透率，在模擬計算時可認為巖塊的滲透系數為0，使得控制方程組中有兩個孔隙度，而只有一個滲透系數。（2）雙重滲透性模型。該模型認為若巖塊中的滲透率不能忽略，則巖塊系統中不僅存在分子彌散作用，而且水流速度不為0。該模型對裂隙和基質分別建立水流方程，主干裂隙建立離散裂隙網絡模型，分枝裂隙連同基質塊建立等效連續介質模型[20]，然后根據邊界通量相同的條件將二者聯系起來。控制方程為：

式中：Ha—主干裂隙水頭/m；

Hr—分枝裂隙水頭/m；

Ka—主干裂隙介質的滲透張量/（m·d?1）；

r—遷移系數；

c—比例系數；

t—時間/d。

顯然，這類模型能較為全面地反映裂隙巖體的滲流特征，比等效連續介質模型更為合理，參數獲取比裂隙網絡模型更為容易，可用來解決區域尺度裂隙水流問題。Barenblatt 等[21]首先提出了“水力雙重介質”的概念，認為裂隙巖石是由空隙性差、導水性強的裂隙系統和空隙性好、導水性差的巖塊孔隙系統共同構成連續介質的統一體，但該模型忽略了裂隙巖體滲流的各向異性特征。Warren 等[22]對裂隙巖體滲透特性的各向異性做了新的假設和描述，但該模型只能應用于均質的正交裂隙網絡。Zimmerman 等[23]針對裂隙/多孔介質的單相流體流動，開發了一種新的雙孔隙度模型，以集中參數方式處理矩陣塊，與Warren-Root 模型對比，節省了近90%的計算時間。楊棟等[24]根據裂隙發育規模與工程尺度的關系，將裂隙巖體看作是由離散介質和擬連續介質組成的廣義雙重介質巖體，提出了廣義雙重介質巖體水力學模型，并對有限元解法進行了較為詳細的研究。雖然該模型在一定程度上刻畫出優先流的現象，但絲毫沒有涉及介質中導水裂隙的具體展布位置，并不能表現出裂隙介質的各向異性、不連續性等特征，模型仍建立在多孔介質滲流理論的基礎上。因此模型應用的缺點有：（1）模型中溶質交換項很難確定，物質交換系數對模型精度影響太大；（2）裂隙網絡不一定能等效為連續介質，適用范圍受限制；（3）求解需迭代計算，較復雜。

1.4 等效-離散耦合模型

該模型是雙重連續介質模型和離散裂隙網絡模型的發展，利用區域分解方法將研究區域按裂隙發育的情況進行分區，不同的分區使用不同的數學模型[25]，即對于裂隙密度大的區域采用等效連續介質模型，對于裂隙密度較小的地區采用離散裂隙網絡模型。該類模型的優點是更符合一般地質條件下裂隙滲流的特征，但也存在水交換量難以確定、兩類模型耦合技術上的問題。Cacas 等[26]運用混合模型對法國某個花崗巖鈾礦不同尺度裂隙中的地下水流過程進行了模擬。黃勇等[27]將基于區域分解法的地下水耦合模型應用于錦屏水電站壩址區三維滲流場的模擬中，該方法有效避免了計算量過大的問題。

2 模型的不確定性分析

由于高放廢物深地質處置深部資料的相對匱乏，加之基巖含水介質高度的非均質和各向異性，進行深地質處置地下水流數值模擬時結果會出現較大的不確定性，因而選用合適的不確定性分析方法進行深地質處置數值模擬不確定性分析顯得尤為重要。

2.1 模型的不確定性來源

地下水流數值模擬的不確定性劃分為三類：

（1）模型參數的不確定性

由于水文地質資料的局限性（包括缺乏資料或參數獲取方法的局限）以及參數在時空上的變異性、尺度效應等都會導致模型參數的不確定性。參數的不確定性分析是影響地下水數值模擬結果可靠性的重要因素。

（2）水文地質概念模型的不確定性

由于資料的缺乏或對水文地質條件認識的局限，建立的水文地質概念模型存在不確定性，包括對地下水流動特征、含水介質的類型、模型邊界及源匯項的估計等。

（3）觀測數據的不確定性

觀測數據的不確定性來源非常廣泛，包括觀測變量隨機分布特征導致的觀測誤差、觀測變量的抽樣誤差、間接測量誤差以及測量儀器本身的觀測誤差和人為記錄誤差等。

以上的不確定性因素，對地下水流模型的穩定性、可靠性和精度有很大影響，因而有必要通過模型的不確定分析，了解模型的不確定因素及其對模型模擬結果的影響程度。

2.2 不確定性分析方法

（1）靈敏度分析

靈敏度分析是一種研究系統各種輸入變化對輸出影響程度的方法。在地下水流數值模擬中，通常以特定的模擬結果為敏感指標（如某個特定點的水位、水質點遷移到特定位置的時間等），以模型的某些要素（如水文地質參數、地下水補排量、邊界條件等）作為敏感因子，改變敏感因子數值，通過模型計算敏感指標的變化程度，分析不同敏感因子的靈敏度。

靈敏度分析方法主要包括全局分析法以及局部分析法。局部分析法是研究某一個敏感因子的變化對敏感指標的影響，該方法主要有：①因子變化法，將敏感因子數值增加或減少一個單位幅度，如5%；②偏差變化法，將敏感因子增加或減少一個標準偏差。全局分析法研究在不同敏感因子的共同作用下，對結果產生的影響，該方法主要有Morris 法、多元回歸法以及Sobol 法等。全局分析法操作相對復雜，目前主要應用于水文模型，在地下水流數值模型靈敏度分析研究中，仍然主要以局部分析法為主。局部靈敏度系數為：

式中：n—分析組數；

yi—參數變化后的輸出值；

y0—參數變化前的輸出值；

Pi—變化后參數值；

P0—初始參數值。

通過靈敏度分析法可將不同敏感因子對地下水流數值模型敏感指標的影響程度進行排序，得到對模擬結果影響較大的敏感因子及其分布區域。此外，通過靈敏度分析獲取更為準確的靈敏度高的敏感因子，可大大提高模型精度，而對于靈敏度不高的敏感因子，則可適當減少獲取該參數的工作量，或者在不確定分析中忽略該敏感因子，達到降維的目的。Mclaren等[28]、Bodvarsson 等[29]在場地尺度上充分考慮裂隙介質的特性，分析了影響地下水運動的參數敏感性，研究了介質對地下水運動的影響和對污染物的阻滯效應。Schwartz[30]利用Toughreact/EOS7R 軟件模擬了裂隙花崗巖中放射性核素遷移過程，就裂隙間距對模擬結果的影響進行了敏感性分析。Shahkarami 等[31]提出一種可以解釋任意長度的衰變鏈以及模擬核素通過任意通道遷移的Sobol 分析方法，并將其應用于模型輸出的全局敏感性分析中。彭志娟[32]分別采用等效連續和雙重介質模型概化裂隙建立高放廢物處置庫開挖擾動區（EDZ）模型，利用Tough2 軟件研究了不同概化模型對129I 和135Cs 在EDZ 中遷移的影響，并對129I 的擴散系數和分配系數進行了敏感性分析。

（2）蒙特卡羅（Monte Carlo）法

蒙特卡羅法是一種被廣泛采用的分析復雜數值模型不確定性分析的方法。該方法是Warren 等[33]最早應用于地下水流研究中，將介質參數作為隨機場，假定已知隨機變量的概率分布和協方差函數，用偽隨機數的生成技術生成成千上萬組輸入變量，產生介質參數隨機分布的大量實現，最終獲得成千上萬組模型計算結果（如地下水模擬中的水頭或溶質濃度）的統計值。該方法回避了隨機分析中的數學困難，只要模擬的次數足夠多，就可得到一個比較精準的概率分布，并且具有收斂速度與問題的維數無關、程序結構簡單等優點。但也存在一定的缺陷，如計算工作量大、誤差難以估計或控制等[34]。因此，目前蒙特卡羅法一般用于比較簡單的模型不確定性分析或者用來驗證其他方法的分析結果。

蒙特卡羅法在國內外的地下水流數值模型不確定性分析中應用十分廣泛。Vaitinen 等[35]利用蒙特卡羅法在芬蘭Romuvaara 場址生成了30 種不同的地下水流模型，用于確定隨機連續介質類型分析的非均勻滲透張量。Joyce 等[36]應用多尺度嵌套建模技術，利用Connectflow 軟件生成了10 個處置庫尺度DFN 模型，分別進行了地下水流動和溶質運移數值模擬。Pandey 等[37]模擬了放射性核素129I 從廢物處理設施釋放并通過地質處置庫遷移到生物圈的過程，利用蒙特卡羅方法進行了阻滯因子、流速以及路徑長度等參數范圍的不確定性分析。蘇銳等[38]以一個假想的CRPGEORC 高放廢物地質處置系統為例，運用Goldsim 軟件對放射性核素79Se 和129I 在該系統遠場中的遷移過程進行了靈敏度分析，同時利用蒙特卡羅法進行了概率統計模擬。林達[39]以某鈾尾礦庫為例，利用GMS軟件結合Monte-Carlo 方法進行了核素U（Ⅵ）在地下水中遷移隨機數值模擬，研究了滲透系數的變異性對尾礦庫地域地下水流動及U（Ⅵ）遷移的影響。

3 數值模擬軟件

目前用于高放廢物深地質處置地下水流數值模擬方面的軟件較多（表1），按照模擬方法的不同，可將軟件分成4 個大類。

表1 主要軟件介紹Table 1 Introduction of main software

（1）等效連續介質模型

Modflow、GMS 以及Porflow 等主要用于孔隙介質和裂隙介質等效連續模型。Belcher[40]利用Modflow-2000 對尤卡山地區進行了比較全面的地下水流預測模擬和分析。L?fman 等[41]通過3D 有限元的模擬方法評估Olkiluoto 場址地下500 m 深部巖體的滲透性，發現施工中的深鉆孔使深部基巖中大量裂隙間的連通性增大。董艷輝等[42]利用Modflow 建立北山地區的地下水流模型，初步劃分了地下水流動系統。王海龍[43]利用GMS 構建了北山區域地下水數值模型，對北山地區地下水流場形態進行了進一步刻畫。季瑞利[44]以北山花崗巖為參考巖性、內蒙古高廟子膨潤土為參考緩沖材料，使用Porflow 對放射性核素135Cs、129I 在處置單元及其圍巖中的遷移行為進行了數值模擬研究。

（2）離散裂隙網絡模型

Connectflow 主要用于離散裂隙網絡的水流、溶質及熱運移模擬。Werner 等[45]采用Connecflow 對瑞典Laxemar 場址高放廢物深地質處置庫建立地下水流數值模型。Joyce 等[36]考慮氣候條件變化，利用Connectflow軟件構建了區域—場址—處置庫多尺度地下水流及溶質運移數值模型。

（3）雙重介質模型

Tough2 系列軟件主要應用于雙重介質的水流、溶質及熱運移模擬。Liu 等[46]、Bodvarsson 等[29]選用Tough2 模擬軟件還原尤卡山在區域尺度上非飽和帶的地下水系統流動規律。Rechard 等[47?48]通過建立溶質運移數值模型驗證尤卡山地區選址的可靠性。Schwartz[30]建立雙孔隙度模型，利用Toughreact 軟件耦合了水流和運移過程，模擬了可能從處置場地中泄露出來的核廢物遷移過程。王禮恒[49]圍繞甘肅北山高放廢物地質處置預選區，分別針對區域—盆地—巖體3 級尺度開展研究，利用Modflow、Tough2 軟件進行3 級尺度地下水流動數值模擬。曹瀟元等[50]采用Tough2-MP/EOS3 和GRACE 重力衛星等方法，建立了北山區域地下水飽和-非飽和流模型。

（4）整合模型

整合模型主要指等效連續介質、離散裂隙網絡、雙重介質以及連續-離散耦合模型的整合，Feflow、Hydrogeosphere 均可用于此類模型的水流、溶質及熱運移模擬。Blessent 等[51]采用Hydrogeosphere 軟件，對區域地下水流場進行了模擬研究，并模擬了不同情景下，裂隙對區域流場的影響以及溶質在裂隙中的運移距離。朱君等[52]應用Feflow 軟件構建三維地下水流數值模型，模擬計算了丘陵山區地下水流動特征下氚的遷移規律。包敏等[53]應用Feflow 建立了熔巖玻璃體239Pu 的溶解釋放和遷移模型，模擬了10 萬年內溶解態239Pu 和膠體態239Pu 的污染羽分布。

4 結論與展望

在高放廢物深地質處置地下水流數值模擬研究中，選取何種合適的模型與對模擬區域的認識程度以及模型需要解決的實際問題有密切的關系。等效連續介質模型適用于大區域、長序列、裂隙發育程度較高或較均勻的地區，方法成熟、所需的數據和參數易于獲取，但是不能精確刻畫裂隙介質中地下水的流動特征。離散裂隙網絡模型適合解決處置場地、儲罐尺度等需要精細刻畫的地下水流問題，但需要大量裂隙及其連通性數據、相關參數等，工作量大、耗時多。雙重介質模型主要用于解決區域尺度裂隙水流問題，但并不能表現出裂隙介質的各向異性、不連續性等特征，適用范圍存在一定的限制。等效—離散耦合模型可以通過區域分解法對裂隙密度大的區域采用等效連續介質模型，對于裂隙密度較小的地區采用離散裂隙網絡模型，更符合一般地質條件下裂隙滲流的特征，但存在交換量難以確定、模型耦合技術問題等。

由于深部構造特征等資料的匱乏、地質屏障的各向異性以及建模要考慮的時間尺度很長，空間尺度差異很大（罐區—處置場—區域水文地質單元），都會導致模擬結果存在較大的不確定性。通過靈敏度分析將不同敏感因子對模型敏感指標的影響程度進行排序，從而提高模型精度、減小參數不確定性分析的工作量。蒙特卡羅法是目前常用的一種模型不確定性分析方法，原理簡單，易于實現。

高放廢物深地質處置工作是一個長期艱巨的任務，目前雖已取得長足的進展，還需注重以下幾方面研究：

（1）數值模型仿真性研究。進一步開展地質、水文地質、裂隙測量等相關的現場調查及監測工作，綜合野外實際觀測資料及地球物理、水文地質手段等獲取裂隙幾何特征及水力學參數，深入開展裂隙發育規律、裂隙滲透性分布特征研究，為裂隙滲流模擬提供準確的裂隙網絡模型和參數，這是提高模型仿真性的基礎。

（2）模型的不確定性分析研究。區域大深度地下水流數值模型存在著較大的不確定性。因而需要根據水文地質結構、參數等建模要素的統計特征和概率分布函數，開展模型的不確定分析，為核廢物處置場選址和安全評價提供服從一定概率分布條件下的要素置信區間。

（3）預測分析模擬研究。預測工況模擬要結合構建場地可能發生的工況，如處置庫長期演變、廢物罐失效情景、極端降雨等情景，開展不同時期、不同場景下的數值模型進行模擬分析。

（4）多介質耦合模型的研究。大區域構建等效連續孔隙介質模型，而場地尺度可以建立裂隙網絡模型，通過兩個模型共同邊界流量和水位相等將兩者耦合，既可精確刻畫場地裂隙流特征，又可宏觀把握水文地質單元地下水流場。