曹瀟元，侯德義，胡立堂

(1.清華大學環(huán)境學院，北京 100084；2.地下水污染控制與修復教育部工程研究中心，北京 100875；3.北京師范大學水科學研究院，北京 100875)

近年來，我國核能事業(yè)取得了飛速發(fā)展，伴隨產(chǎn)生了大量放射性廢物，包括高、中、低放廢物，其中高放廢物具有衰變時間長、釋熱率高、放射性強的特點，對人類生存環(huán)境與生命安全存在潛在威脅。因此，高放廢物能否得到安全的處理與處置，將制約著我國核工業(yè)的發(fā)展，也是人類安全與環(huán)保方面亟待解決的重大課題[1]。2017年，我國發(fā)布的《中華人民共和國核安全法》明確規(guī)定高水平放射性廢物應實行集中深地質(zhì)處置。在高放廢物的深地質(zhì)處置方面，美國、瑞典、芬蘭、法國等目前在該應用領域和研究處于世界先進水平，同時也為我國高放廢物處置的發(fā)展提供了重要參考[2]。Montazer等人建立了尤卡山非飽和帶地下水模型，為高放核廢料選址提供了一個初步的模型框架[3]。隨后，一些專家學者對尤卡山地區(qū)的非飽和帶-飽和帶中的地下水流動以及溶質(zhì)運移開展了大量數(shù)值模擬研究[4-5]，其中包括尤卡山地區(qū)非飽和帶的地下水流研究以及多孔介質(zhì)中溶質(zhì)遷移及裂隙流等研究。我國高放核廢料處置的研究起步較晚。1985—1986年，經(jīng)過全國地質(zhì)資料的綜合分析對比，從地質(zhì)條件、地質(zhì)構造、巖性以及水文地質(zhì)條件、自然地理、經(jīng)濟地理等因素考慮，在全國范圍內(nèi)篩選出五個候選區(qū)：華東預選區(qū)、華南預選區(qū)、西南預選區(qū)、內(nèi)蒙預選區(qū)和甘肅北山預選區(qū)，處置庫圍巖包括花崗巖、凝灰?guī)r、泥巖和頁巖。至1989年，我國高放廢物處置庫選址工作的重點主要集中在甘肅北山地區(qū)。

掌握區(qū)域地下水動力特征是高放廢物地質(zhì)處置場選址及適宜性評價的基礎條件。因此，建立區(qū)域地下水流模型是高放廢物地質(zhì)處置研究中最基礎和最關鍵的工作之一，它涵蓋整個高放廢物地質(zhì)處置區(qū)，具有宏觀掌握區(qū)域地下水流場形態(tài)的作用，模型的模擬結果可以為預選區(qū)提供合理的概念模型，并為處置庫的選擇、評價、材料選取、處置方案設計等問題提供基礎條件，也是論證候選場址可行性和對各種技術方案有效預測的重要依據(jù)[6-7]。國內(nèi)一些專家學者對甘肅北山預選區(qū)進行地質(zhì)試驗[8-11]，主要對北山地區(qū)地下水流動介質(zhì)的滲透性以及研究區(qū)地下水的循環(huán)屬性進行研究。隨后一些學者開展了模擬研究[12-15]。但是，由于區(qū)域地下水流運動模擬涉及的尺度大、參數(shù)多，傳統(tǒng)地下水數(shù)值模擬軟件(如MODFLOW、GMS等)不能反應多組分和多相流的地下水流動特征，并且缺乏區(qū)域尺度上對地下水流的精細刻畫。

本研究為了保證模型計算的精度與可靠性，利用國際上多相多組分流體的并行版數(shù)值模擬軟件TOUGH2-MP，靈活、精細地刻畫了干旱地區(qū)區(qū)域地下水的補給類型以及地下水流場的分布形態(tài)，基于已有的水文地質(zhì)調(diào)查成果，建立數(shù)值模型，分析模擬與實測地下水位數(shù)據(jù)的擬合情況，將地下水總儲量的變化趨勢作為模型的識別項，驗證模型參數(shù)的合理性以及準確性，為精細劃分區(qū)域地下水流動系統(tǒng)提供可靠的技術支撐，為預選區(qū)選址的安全性評價提供重要依據(jù)。

1 研究區(qū)概況

研究區(qū)位于甘肅省河西走廊以北，南距玉門鎮(zhèn)60 km，北部延伸至與蒙古的國界，東部包括黑河流域(圖1)，研究總面積約69 125 km2。在區(qū)域地勢上具有西高東低、北高南低的特點，地形起伏較大。研究區(qū)標高在西部約2 300 m，由北至南的河西走廊一帶地面標高逐漸降低到約1 200 m。從西部的馬鬃山至東部黑河流域，研究區(qū)的地面標高降低到約900 m，成為最低區(qū)。研究區(qū)分布有大面積侵入花崗巖，上覆少量第四系沉積物以及黏土夾層。由于花崗巖受多期構造運動影響，發(fā)育有斷層和裂隙，主要以東西走向的壓性逆斷層為主。研究區(qū)為我國典型的干旱氣候區(qū)，具有降水量少、蒸發(fā)量大、干冷多風的特點，多年平均降水量為60～80 mm，降水量的60%以上集中在6—8月份，年平均蒸發(fā)量3 200 mm，年平均氣溫4～7 ℃，年平均風速達4 m/s。

通過已獲得的資料和數(shù)據(jù)分析，研究區(qū)含水介質(zhì)可分為裂隙介質(zhì)和孔隙介質(zhì)兩種類型。地下水系統(tǒng)可分為三個子系統(tǒng)，即基巖裂隙水子系統(tǒng)、孔隙-裂隙水子系統(tǒng)、孔隙水子系統(tǒng)三大類。此外，根據(jù)中國1∶250萬地質(zhì)圖推斷可得，北山地區(qū)涵蓋了3條明顯斷層(圖1)。

圖1 研究區(qū)示意圖Fig.1 Schematic map of the study area

2 研究方法

TOUGH2是目前國際上可以模擬多組分、多相流的數(shù)值模擬程序之一，相比其他同類模擬軟件，它可以較靈活處理非均勻介質(zhì)體及裂隙等問題的模擬。TOUGH2-MP的研發(fā)實現(xiàn)了大規(guī)模飽和-非飽和帶地下水流模型的精細計算。TOUGH2程序采用積分有限差分法對網(wǎng)格進行空間離散化求解，僅需要通過單元及相鄰單元間的局部幾何信息，即可實現(xiàn)對網(wǎng)格的任意剖分。對于時間離散化，程序對任意時間步長采用牛頓-拉弗森方法對質(zhì)能守恒方程的殘差方程組進行反復迭代求解，直到達到收斂標準。

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

模型頂界面為大氣壓邊界，取為恒溫(15 ℃)恒壓(標準大氣壓)邊界。模型的底面概化為隔水邊界。根據(jù)美國地質(zhì)調(diào)查局NASA的航天飛機雷達地形測繪數(shù)據(jù)(SRTM)，獲得研究區(qū)的高程，結合高程及河網(wǎng)分布的情況，將研究區(qū)邊界分為流量邊界和隔水邊界。由于在TOUGH2軟件中對流量邊界難以給定，因此在模型中將流量邊界定為定壓邊界(圖1中的Q1和Q2邊界)。

根據(jù)地質(zhì)大調(diào)查以及核工業(yè)北京地質(zhì)研究院的實時補充調(diào)查成果[16]，研究區(qū)主要巖性是變質(zhì)巖和花崗巖，部分區(qū)域上伏有約2～15 m的第四系沉積物，含水層介質(zhì)包括孔隙和裂隙。變質(zhì)巖及侵入巖的滲透系數(shù)非常低，然而，研究區(qū)分布有裂隙和斷層，因此，在裂隙和斷層的交匯處，滲透系數(shù)相對較高。研究區(qū)的變質(zhì)巖及侵入巖可以分為3類：強風化巖層、中風化巖層和弱風化巖層，厚度分別為50，100，500 m。為了更清楚地展現(xiàn)地下水位變化以及非飽和帶的地下水飽和度變化，模型中，將強風化巖層的50 m細化為15層，厚度分別為1.0，1.4，0.5，1.0，1.0，2.0，2.0，2.0，2.5，2.5，3.0，3.0，9.0，10.0，10.0 m。中等風化巖層可以細化為5個等厚巖層，厚度為20 m。弱風化巖層分為3層，分別是100，200，200 m。模型總共分為23層。

2.2 模型控制方程

EOS3模塊可以應用于非飽和帶—飽和帶區(qū)域地下水流和熱運移的模擬計算，該模塊已經(jīng)成功地應用在許多高放核廢料處置庫選址的工程中[17-18]，該模塊的質(zhì)量與能量控制方程為：

質(zhì)能控制方程：

質(zhì)量累計量：

質(zhì)量對流通量：

能量對流通量：

F3=-λ

能量累積量：

相速：

式中：Mκ——區(qū)域內(nèi)每單位體積的質(zhì)能累計量/(kg·m-3)；

κ——成分,取值1(水)，2(氣)，3(溫度)；

Vn——第n個網(wǎng)格單元的體積/m3；

F——力/Pa；

n——內(nèi)部的單位法線向量；

Γn——封閉曲面的面積/m2；

qκ——區(qū)域內(nèi)每單位體積的質(zhì)源或能源/(kg ·m-3·s-1)；

φ——孔隙介質(zhì)的孔隙率；

Sβ——各個相的飽和度；

ρβ——各個相的密度/(kg·m-3)；

β——相態(tài)，取值G(氣相)，L(液相)；

uβ——各個相中的達西流/(m·s-1)；

λ——導熱率/(W· K-1· m-1)；

μβ——黏度/(Pa·s)；

hβ——在各個相中的比熱焓/(J·kg-1)；

CR——導熱系數(shù)/(W ·K-1·m-1)；

T——溫度/(℃或者K)；

krβ——相對相態(tài)的滲透率；

g——重力加速度/(m·s-2)；

k——絕對滲透率/m2。

2.3 數(shù)值模型

2.3.1高程及網(wǎng)格離散高程對地下水流場的分布情況起到非常重要的作用。在本研究中，利用美國NASA宇航局的衛(wèi)星數(shù)據(jù)(SRTM)得到了高程數(shù)據(jù)。研究區(qū)的高程西高東低，最高點在西部的馬鬃山，海拔約2 400 m。研究區(qū)東北角黑河流域海拔較低，最低點可達400 m。本次研究利用胡立堂等[19]開發(fā)的IGMESH軟件。該軟件是一款專門針對TOUGH2編制的前后處理軟件，可用于windows 98/2000/XP/Vista/windows 7等微軟主流的操作系統(tǒng)，具有可視化操作界面。利用IGMESH軟件將研究區(qū)在平面上劃分為12 390個不規(guī)則的網(wǎng)格單元?？偩W(wǎng)格數(shù)為284 970個，總鏈接數(shù)為1 120 107個(圖2)。為了更精細地刻畫高放廢物處置的預選區(qū)的地下水流場，在預選區(qū)位置對網(wǎng)格進行剖分。

圖2 研究區(qū)網(wǎng)格剖分及高程Fig.2 Mesh and elevation of the study area

2.3.2主要水文地質(zhì)參數(shù)

為探究研究區(qū)滲透系數(shù)隨著深度增加的變化情況，參照核工業(yè)北京地質(zhì)研究院的野外調(diào)查數(shù)據(jù)[16]，分析研究區(qū)水文地質(zhì)參數(shù)的數(shù)據(jù)分別有：①我國第一輪水文地質(zhì)普查的37個抽水試驗點；②本研究團隊通過雙環(huán)入滲試驗獲得淺層井121個試驗數(shù)據(jù)；③核工業(yè)北京地質(zhì)研究院采用雙栓塞水文地質(zhì)試驗系統(tǒng)獲得的兩口深鉆孔BS05與BS06。其中，通過抽水試驗，測得淺層井的滲透系數(shù)變化范圍0.001 5～3.82 m/d。利用雙環(huán)入滲試驗獲得的淺層121組試驗數(shù)據(jù)為0.11～18.8 m/d，平均滲透系數(shù)約4.27 m/d。深鉆孔BS05、BS06的有效試驗深度為44.22～600.55 m，BS05與BS06的平均滲透系數(shù)分別為0.008 8 m/d、0.000 18 m/d。通過數(shù)據(jù)結果整理分析可知，研究區(qū)內(nèi)風化巖的滲透系數(shù)隨著深度增加而衰減，模型中的平均滲透率為1×10-13～1×10-17m2。在本次研究中，詳細的水文地質(zhì)參數(shù)分配見表2。在模型中，巖石顆粒密度、巖層熱傳導率、孔隙壓縮比分別為2 650 kg/m3，1 W/(m·℃)，3.7×10-10Pa-1。頂層的巖石比熱設置為8.0×1099J/(kg·℃)，其他層位設置為1 000 J/(kg·℃)。

研究區(qū)存在區(qū)域斷層，主要以東西向延展為主，深度可延伸至約80 m。斷層中的裂隙被花崗巖與石英填充，經(jīng)2013—2014年的裂隙統(tǒng)計數(shù)據(jù)分析，多數(shù)裂隙的走向主要集中分布在NW(300°～330°)以及NE(40°～70°)，絕大多數(shù)的裂隙傾角集中在60°～70°，也有小部分裂隙傾角在40°～50°，90%以上的地表裂隙開度大于10 mm，長度多超過1 m。在水文地質(zhì)參數(shù)分區(qū)中，將區(qū)域裂隙單獨分區(qū)(表1)，孔隙度設為0.30。

表1 研究區(qū)水文地質(zhì)參數(shù)分區(qū)和參數(shù)值Table 1 Hydrogeology zones and parameters for the 23 model layers

注：分區(qū)中2-9分別代表第四系沉積物、黏土夾亞砂土、砂土以及不同風化程度的花崗巖；10代表裂隙。kx，ky，kz分別表示x、y、z方向的滲透率。

2.3.3初始條件及源匯項

模型的初始流場是在降水90 mm條件下，模型長時間運算后達到靜水壓力平衡狀態(tài)下的初始狀態(tài)。在本次研究中，地表溫度假設在15 ℃，并按照地溫梯度30 ℃/km計算溫度分布。本次模型計算中，不考慮溫度隨時間的變化。水頭與壓力的關系為：

(1)

式中：p——地下水的壓強/Pa；

p0——標準大氣壓/Pa；

ρ——水的密度/(kg·m-3)；

g——重力加速度/(m·s-2)；

H——水頭/m；

z——位置高度/m。

模型地下水補給來源包括溝谷洪流入滲和降水入滲，地下水的排泄方式包括蒸發(fā)排泄、泉水出露及以側向徑流的方式向區(qū)外排泄，在東部部分區(qū)域排入黑河沖積平原。由于北山地區(qū)沒有氣象站，本次研究借鑒鄰近研究區(qū)和全球降水氣候計劃(Global Precipitation Climatology Project，簡稱GPCP)衛(wèi)星的降水數(shù)據(jù)，得到2003—2014年的月降水量數(shù)據(jù)。從數(shù)據(jù)中可以看出，北山地區(qū)降水量很小，6—9月是降水量最大的時期，最大月降水量可達47 mm，月降水量最小只有0.89 mm。研究區(qū)多年平均蒸發(fā)量從西向東分為3個帶：東部蒸發(fā)強度為3 500 mm/a，中部蒸發(fā)強度為3 100 mm/a，西部蒸發(fā)強度為2 800 mm/a，多年平均蒸發(fā)強度為2 900～3 200 mm/a。

3 模型識別

本次研究的模型識別期間為2003—2014年，分別利用鉆孔地下水位對比、地下水的總儲量變化趨勢作為模型的識別手段，證明模型參數(shù)符合實際區(qū)域的水文地質(zhì)條件，并可以比較準確地反應出區(qū)域地下水流動特征、分布形態(tài)以及巖性的滲流特征，驗證了模型的準確性以及可靠性。

3.1 統(tǒng)測和模擬水位對比

由于研究區(qū)內(nèi)沒有長時間序列的觀測值，僅收集到2008—2009年33口統(tǒng)測井的地下水位。因此，將已有實測水位與計算水位進行對比。約69 125 km2的模擬面積上91%的模擬和觀測的地下水位絕對誤差在5 m以內(nèi)，絕對誤差在1 m以內(nèi)的占64%，擬合結果較為合理，模擬結果可以反映區(qū)域地下水流的特征(圖3)。

圖3 計算與實測地下水位數(shù)據(jù)對比圖Fig.3 Fitting between the observed data and the simulated data

3.2 地下水儲量變化與GRACE衛(wèi)星反演數(shù)據(jù)對比

由于研究區(qū)面積較廣且鉆孔的位置相對集中，因此，僅通過擬合鉆孔的實測值與模擬值并不能完全反映模型的可靠性與準確性。GRACE衛(wèi)星可以通過重力場的變化反演推算地下水的變化。此外，由于該數(shù)據(jù)具有成本低、更新速度快的特點，因此GRACE衛(wèi)星數(shù)據(jù)在國際上已取得較廣泛的應用[12-14]。在本次研究中，利用GRACE重力衛(wèi)星和GLDAS陸面水文模型反映出該區(qū)域的地下水儲量變異，然后將其相對變化與模擬的地下水儲量變化進行比較。圖4反映了2003—2009年GRACE反演的地下水儲量變異量與計算的地下水儲量變化的對比情況。模型計算結果可以反映區(qū)域地下水儲量呈弱下降趨勢。模型計算結果與GRACE衛(wèi)星反演數(shù)據(jù)的變化趨勢一致。

圖4 GRACE衛(wèi)星反演的陸地水儲量數(shù)據(jù)與模型計算值對比Fig.4 Comparison between the simulated and GRACE-derived TWS variations

4 結果分析

4.1 地下水均衡分析

降雨是研究區(qū)地下水的主要補給來源。由前期調(diào)查和相關研究，研究區(qū)入滲系數(shù)小，大部分降雨通過地表徑流形成溝流，經(jīng)過蒸發(fā)或地面徑流流出區(qū)外。在降水比較集中的夏秋季節(jié)，部分降水以滲流的方式，沿著構造帶或者第四系松散沉積物出露地帶向下補給地下水。但由于研究區(qū)內(nèi)大部分區(qū)域花崗巖出露，因此，比較集中的降水多以地表徑流的方式沿著溝谷洼地向地勢較低的地方排泄，并在地表形成大量的沖溝。為了研究降水補給地下水部分和地表徑流排泄部分各自所占比例，研究區(qū)水資源的均衡計算：

A×P=ΔW+Q1+Q2+QS

(2)

式中：P——某年的降水量/m；

A——研究區(qū)面積/m2；

ΔW——研究區(qū)某年地下水儲量變化值/m3，正值表示增加，負值表示衰減；

Q1、Q2——東流量邊界及南流量邊界的流出量/m3，根據(jù)達西公式由模型邊界網(wǎng)格計算得出；

QS——某年的地表徑流排泄量和蒸發(fā)量之和/m3。

研究區(qū)地下水儲量變化可由TOUGH2模型計算得出。研究區(qū)蒸發(fā)及地表徑流量通過式(2)進行間接計算獲得。為研究極端氣候條件對研究區(qū)地下水的影響，本研究考慮了不同極端氣候條件(200 mm和300 mm降水情景)模擬區(qū)的地下水資源均衡狀況，結果見表3。從表中可知，現(xiàn)狀(90 mm降水情景)條件下，大氣降雨的補給量約為0.655×108m3，地下水側向年排泄量約為103.12×104m3(Q1段為30.76 kg/s；Q2段為1.94 kg/s)，模擬的地下水儲量減少量約為32.2×104m3，占降水補給量的0.49%，該比例可能是模型計算誤差引起，說明整個研究區(qū)地下水資源處于均衡狀態(tài)。在極端降雨事件情景(200，300 mm)時，模擬一萬年的地下水側向排泄量和地下水儲量減少量變化并不大，地下水資源也處于均衡狀態(tài)(計算的儲量變化分別占降水補給量的0.23%和0.15%)。結果表明，99%左右的降水量通過地表徑流和蒸發(fā)被排泄，研究區(qū)只有1%左右的降水量入滲至地下(表2)。

表2 不同降雨情景模擬的地下水均衡情況Table 2 Variation in groundwater budgets in three scenarios

4.2 區(qū)域地下水系統(tǒng)劃分

相比地表水流動系統(tǒng)，地下水流動系統(tǒng)并沒有明顯的分界線，然而科學地劃分區(qū)域地下水流動系統(tǒng)是開展地下水相關研究的基礎，也是正確認識地形地貌、地層巖性以及構造控制條件下地下水演化的必然結果，對分析研究區(qū)的地下水流動特征、高放廢物地質(zhì)處置選址與安全性評價都有著重要的指導作用。

圖5展示了北山地區(qū)地下水流場的分布情況。研究區(qū)中的地下水系統(tǒng)由北向南可依次分為三個獨立的水文地質(zhì)單元(Ⅰ、Ⅱ、Ⅲ)。不同的地下水循環(huán)系統(tǒng)有不同的地下水流向及循環(huán)特征。水文地質(zhì)單元Ⅰ位于研究區(qū)北部，面積約為研究區(qū)面積的1/2，南部邊界主要集中在研究區(qū)中部的山脊一帶。受地勢因素影響，水文地質(zhì)單元Ⅰ的地下水在西部匯集最終流向東部邊界。水文地質(zhì)單元Ⅱ位于研究區(qū)南部，面積約為研究區(qū)面積的1/3,北部與水文地質(zhì)單元Ⅰ相鄰，南部邊界主要集中在變質(zhì)巖與黏土及砂壤土的交界處。該水文地質(zhì)單元與水文地質(zhì)單元Ⅰ的地下水流向基本一致，兩者最終匯集到黑河流域下游平原區(qū)。水文地質(zhì)單元Ⅲ位于研究區(qū)最南部，該水文地質(zhì)單元的巖性多以變質(zhì)巖為主，且地下水向南邊界流動，最終排入河西走廊地區(qū)。

高放廢物處置庫的預選區(qū)完全處于水文單元Ⅲ中，由此可知，通過本次模擬結果可以很大程度上縮小研究區(qū)范圍，并可根據(jù)驗證后的水文地質(zhì)單元參數(shù)，對預選區(qū)的核素遷移模擬進行應用。在該研究區(qū)中，可以清晰地發(fā)現(xiàn)研究區(qū)的地下水位分布與自然條件下河網(wǎng)的分布形態(tài)基本保持一致，符合西高東低的地形地貌條件，并與地表水流動方向保持一致，該模擬結果與肖豐提出的北山地區(qū)地下水徑流特征與溝系分布形態(tài)總體一致的結論相符[10]。一方面證明了該模型模擬結果的可靠性，另一方面表明研究區(qū)地形是決定地下水流動的主要影響因素。

圖5 北山地區(qū)地下水流場分布(2014年12月)Fig.5 Groundwater flow field in the Beishan area at the end of 2014

5 結論

我國甘肅北山地區(qū)屬于極干旱的戈壁地區(qū)，在有限的經(jīng)濟投入的前提下，了解區(qū)域地下水動力特征是地下工程建設的前提條件，而數(shù)值模擬是推斷研究區(qū)水文地質(zhì)條件的重要手段。本次研究利用國際上先進的多相流數(shù)值模擬軟件TOUGH2-MP/EOS3模塊對研究區(qū)的地下水流動進行模擬分析，并根據(jù)區(qū)域的地形地貌、地下水循環(huán)交替特征，以及介質(zhì)的巖性特征等因素，對研究區(qū)的地下水流動系統(tǒng)進行了劃分，該方法在針對數(shù)據(jù)資料少、研究區(qū)面積廣的區(qū)域地下水流場模擬刻畫方面具有較好的適用性及有效性。獲得的主要結論為：

(1)利用數(shù)值模擬方法，探討了區(qū)域地下水動力特征不明確、流場難以準確描述的問題。通過地下水均衡分析發(fā)現(xiàn)，研究區(qū)約99%的降水量通過地表徑流和蒸發(fā)被排泄掉，1%的降水量入滲至地下。

(2)通過調(diào)整模擬過程中的不同氣候條件，發(fā)現(xiàn)由于研究區(qū)的巖性具有較低的入滲率，導致不同氣候條件對區(qū)域地下水流場的影響較小；根據(jù)區(qū)域地下水系統(tǒng)劃分可知，研究區(qū)地形是影響地下水流動的主要影響因素。

(3)根據(jù)地下水流場的模擬形態(tài)，可將區(qū)域地下水流動系統(tǒng)分為三個子系統(tǒng)(Ⅰ、Ⅱ、Ⅲ)，且預選區(qū)正位于區(qū)域內(nèi)獨立的地下水子系統(tǒng)Ⅲ中。該研究成果可為進一步分析核素在預選區(qū)的遷移范圍提供基礎。