張 雪，胡繼華，施曉文，王小元

(中國(guó)核電工程有限公司，北京 100840)

為保證核設(shè)施建設(shè)、運(yùn)行及退役等階段土壤及地下水環(huán)境的安全，根據(jù)國(guó)家核安全法規(guī)、標(biāo)準(zhǔn)和導(dǎo)則的要求，必須通過(guò)實(shí)驗(yàn)和數(shù)值模擬等手段，充分掌握核設(shè)施廠址地區(qū)放射性核素在土壤及水環(huán)境中的遷移擴(kuò)散規(guī)律，提供工程設(shè)計(jì)所需的參數(shù)資料。

以往對(duì)假想事故條件下放射性液體的地下水環(huán)境影響評(píng)價(jià)，假設(shè)放射性液體瞬時(shí)進(jìn)入飽和地下水環(huán)境，不考慮包氣帶對(duì)核素遷移的遲滯作用。對(duì)于放射性核素來(lái)說(shuō)，包氣帶對(duì)放射性核素的阻滯作用非常明顯，李書(shū)紳等[1]的研究表明，在自然條件下，可吸附核素的遷移距離非常短，因此研究放射性核素在包氣帶中遷移是必要的。放射性核素在包氣帶中的遷移主要以垂直向?yàn)橹鳎遗c地下水環(huán)評(píng)的評(píng)價(jià)范圍相比，水平向的運(yùn)動(dòng)幾乎可以忽略。

將包氣帶水分運(yùn)動(dòng)處理為一維垂向運(yùn)動(dòng)，飽和帶地下水為三維運(yùn)動(dòng)，將兩種不同維數(shù)運(yùn)動(dòng)的問(wèn)題進(jìn)行耦合而得到擬三維數(shù)值模型。包氣帶一維垂向模擬選用Hydrus-1D軟件，飽和帶三維模擬選用GMS軟件。將潛水面作為包氣帶模型下邊界，同時(shí)作為飽和帶模型上邊界，通過(guò)包氣帶模擬獲得包氣帶下邊界地下水流量及放射性核素濃度隨時(shí)間變化值，并將其作為飽和帶的源匯項(xiàng)(上邊界條件)輸入到飽和帶模型，實(shí)現(xiàn)包氣帶-飽和帶的松散耦合。

1 國(guó)內(nèi)外研究進(jìn)展

大多數(shù)飽和-非飽和水流耦合模型都基于MODFLOW建立，基于一維非飽和水流運(yùn)動(dòng)和三維地下水運(yùn)動(dòng)的含水層數(shù)值模擬系統(tǒng)如圖1所示。

圖1 基于MODFLOW建立的擬三維飽和-非飽和水流系統(tǒng)

Havard等人[2]建立了LINKFLOW模型，該模型將模擬非飽和區(qū)垂向流動(dòng)的Richards方程的有限差分解與模擬飽和區(qū)水流的三維有限差分模型MODFLOW進(jìn)行耦合。Facchi等人(2004)[3]將土壤水運(yùn)動(dòng)模型SVAT與MODFLOW進(jìn)行了耦合。Niswonger和Prudic (2004)[4]、Niswonger等人(2006)[5]在MODFLOW-2005 的基礎(chǔ)上發(fā)展了非飽和水流運(yùn)動(dòng)計(jì)算程序包UZF1。Twarakavi等人(2008)[6]采用類(lèi)似的方法建立了一個(gè)更為完善的飽和-非飽和水流運(yùn)動(dòng)模型，采用HYDRUS-1D程序包模擬非飽和帶水流運(yùn)動(dòng)，并與MODFLOW-2000進(jìn)行耦合。Yakirevich等人(1998)[7]采用有限差分法建立了擬三維飽和-非飽和水流運(yùn)動(dòng)和溶質(zhì)運(yùn)移耦合模型。Kuznetsov等人(2012)[8]進(jìn)一步發(fā)展了Yakirevich的耦合模型，采用完全三維水流運(yùn)動(dòng)方程描述地下水與毛細(xì)管區(qū)水流運(yùn)動(dòng),并與非飽和一維水流運(yùn)動(dòng)方程進(jìn)行耦合。

林琳等(2005)[9]根據(jù)多孔介質(zhì)中區(qū)域飽和-非飽和水分運(yùn)動(dòng)的特征，建立基于迦遼金有限單元法的擬三維模型。朱炎(2013)[10]將非飽和帶水流運(yùn)動(dòng)簡(jiǎn)化為一維垂向運(yùn)動(dòng)，保持飽和帶三維運(yùn)動(dòng)的特點(diǎn)，將兩者進(jìn)行耦合。查元源(2014)[11]將Richards方程模型及水均衡模型與三維地下水模型耦合，建立了區(qū)域飽和-非飽和三維水分運(yùn)動(dòng)耦合模型。于洋等(2017)[12]通過(guò)GMS和Hydrus的結(jié)合，實(shí)現(xiàn)了包氣帶-飽和帶水鹽運(yùn)移的耦合模擬。張雅楠(2019)[13]明晰了HYDRUS-MODFLOW模型的耦合原理及運(yùn)行機(jī)制。

2 研究區(qū)概況

研究區(qū)屬大陸性干旱氣候，降水少，蒸發(fā)強(qiáng)烈。地勢(shì)北高南低，北部為基巖山地，海拔1 250.0～1 440.0 m，南部為碎石戈壁，海拔一般為1 229.0～1 275.0 m。多年平均降雨量62.9 mm，多年平均蒸發(fā)量為2 338.9 mm。地下水可分為松散巖類(lèi)孔隙水和基巖裂隙水兩類(lèi)。松散巖類(lèi)孔隙水區(qū)域滲透系數(shù)為2.78×10-4～7.20×10-3cm/s；基巖裂隙水區(qū)域滲透系數(shù)為3.97×10-6～1.53×10-4cm/s。地下水水位埋深在2.80 m～20.06 m之間，水位標(biāo)高在1 205.25 m～1 227.8 m之間。

3 包氣帶模擬研究

3.1 數(shù)值模型的建立

包氣帶一維垂向模擬選用Hydrus-1D軟件，建立模型時(shí)做如下假設(shè)：包氣帶中水流和核素為垂向一維流，核素的釋放源在地表；包氣帶介質(zhì)對(duì)核素的吸附是可逆的，為一級(jí)線性反應(yīng)；只存在一級(jí)化學(xué)反應(yīng)，且不考慮核素的衰變產(chǎn)物；忽略溫度和氣體的影響。

試驗(yàn)場(chǎng)地包氣帶厚度約3 m，假設(shè)罐體中的放射性液體發(fā)生泄漏，并由地表垂直入滲補(bǔ)給包氣帶，泄漏總量為8.50 m3，核素種類(lèi)包括C-14、Se-79、Pd-107、Cs-135和Sn-126五種，活度濃度分別為9.40×104Bq/L、8.80×103Bq/L、1.50×104Bq/L、5.60×104Bq/L和6.64×104Bq/L。

3.2 初始條件

(1)假設(shè)初始條件下水流處于穩(wěn)定狀態(tài)，水流模型初始值由壓力水頭給出：h(z,0)=h0(z)，-H≤z≤0。式中：h為土壤水壓力冰頭，H為包氣帶厚度，z為垂直坐標(biāo)，設(shè)潛水面處壓力水頭為0，向上遞減，地表處壓力水頭為-3 m。

(2)核素遷移模型初始時(shí)刻選在放射性液體剛從地表泄漏，包氣帶中放射性液體初始濃度設(shè)為0。

3.3 邊界條件

(1)水流模型上邊界為變水頭/流量邊界:

(1)

泄漏會(huì)產(chǎn)生地表積水，t1為泄漏持續(xù)時(shí)間，H(t)為地表積水厚度。即地表有放射性液體泄漏時(shí)(0≤t≤t1)，上邊界為隨時(shí)間變化的水頭邊界；泄漏停止后(t>t1)上邊界為0流量邊界。k(h)為非飽和水力傳導(dǎo)度。

下邊界選擇潛水面，設(shè)為定水頭邊界：

h(z,t)=0，z=-H，t≥0

(2)

(2)核素遷移模型上邊界為定濃度通量邊界：

(3)

式中，θ為體積含水率，D為水動(dòng)力彌散系數(shù)，v為水流垂向流速，c為核素在包氣帶中的濃度，qs(t)為地表水分通量，C為泄漏液體中放射性核素濃度。

下邊界為零濃度梯度邊界。

3.4 模型主要參數(shù)確定

土-水特征曲線參數(shù)θs、θr、n、α和Ks通過(guò)現(xiàn)場(chǎng)土-水特征曲線試驗(yàn)和現(xiàn)場(chǎng)包氣帶滲水試驗(yàn)獲取初值，并利用包氣帶水分遷移試驗(yàn)，通過(guò)模型識(shí)別驗(yàn)證獲取模型參數(shù)取值；干密度ρd取測(cè)試值；彌散度aL選擇多孔示蹤試驗(yàn)值。參數(shù)取值列于表1。

表1 參數(shù)取值表

表2 分配系數(shù)和衰變常數(shù)取值

3.5 模擬結(jié)果及分析

根據(jù)包氣帶-飽和帶耦合模擬思路，包氣帶的模擬結(jié)果將作為飽和帶的源匯項(xiàng)(上邊界條件)輸入到飽和帶模型中，模擬得到潛水面處水流速度及C-14、Se-79、Pd-107、Cs-135和Sn-126五種核素濃度隨時(shí)間變化如圖2～圖7所示。

圖2 潛水面處水流速度隨時(shí)間變化

圖3 潛水面處C-14活度濃度隨時(shí)間變化

圖4 潛水面處Se-79活度濃度隨時(shí)間變化

圖5 潛水面處Pd-107活度濃度隨時(shí)間變化

圖6 潛水面處Cs-135活度濃度隨時(shí)間變化

圖7 潛水面處Sn-126活度濃度隨時(shí)間變化

由圖2～7模擬數(shù)據(jù)結(jié)果可知：

(1)包氣帶對(duì)水流運(yùn)動(dòng)(核素遷移)存在延遲作用：放射性液體并非瞬時(shí)進(jìn)入飽和帶，包氣帶具有一定的滲透能力，水流只能以一定的速率滲入飽和地下水環(huán)境(本次模擬完成水流入滲大約需要30 d)，核素主要在對(duì)流作用的驅(qū)動(dòng)下滲入飽和地下水環(huán)境，也即包氣帶對(duì)核素遷移存在延遲作用。

(2)包氣帶對(duì)核素遷移存在滯留作用：由于包氣帶介質(zhì)對(duì)核素的吸附作用，使遷移至潛水面處核素的活度濃度緩慢增大而并非瞬時(shí)增大。對(duì)于吸附性較強(qiáng)的核素(Se-79、Pd-107、Cs-135和Sn-126)，其遷移至潛水面處的核素活度濃度遠(yuǎn)小于泄漏的放射性液體核素活度濃度，Se-79、Pd-107、Cs-135和Sn-126在潛水面處的活度濃度峰值分別為3.67 Bq/L、705 Bq/L、3.67 Bq/L和0.045 Bq/L，對(duì)應(yīng)的年待積有效劑量分別為8.19×10-4mSv、0.018 9 mSv、8.19×10-4mSv和5.83×10-5mSv，已小于WHO規(guī)定的0.1 mSv。而對(duì)于吸附性較弱的C-14，遷移至潛水面處的核素活度濃度與泄漏液核素活度濃度差別不大，C-14在潛水面處的峰值活度濃度為9.37×104Bq/L，為泄漏液活度濃度的99.68%，對(duì)應(yīng)的年待積有效劑量為15.50 mSv，與泄漏的放射性液體(15.55 mSv)相比只降低了0.05 mSv。

即包氣帶對(duì)核素的遷移存在明顯的遲滯作用，根據(jù)包氣帶模擬結(jié)果，接下來(lái)只對(duì)潛水面處核素活度濃度對(duì)應(yīng)年待積有效劑量大于評(píng)價(jià)準(zhǔn)則(0.1 mSv)的C-14開(kāi)展耦合模擬研究。

4 飽和帶地下水水流數(shù)值模型

4.1 模型范圍、邊界條件與水力特征

模擬范圍與邊界條件充分考慮評(píng)價(jià)區(qū)水文地質(zhì)條件確定。

模擬區(qū)內(nèi)地下水運(yùn)動(dòng)符合達(dá)西定律。地下水系統(tǒng)的垂向運(yùn)動(dòng)主要是大氣降水在地表入滲后向下補(bǔ)給含水層，孔隙水含水層垂向補(bǔ)給基巖裂隙水，模擬區(qū)內(nèi)沒(méi)有明顯的隔水層，根據(jù)地下水位長(zhǎng)期監(jiān)測(cè)結(jié)果，該地區(qū)地下水水位隨季節(jié)變化不明顯。因此，模擬區(qū)地下水系統(tǒng)可概化為非均質(zhì)三維穩(wěn)態(tài)地下水流系統(tǒng)。

4.2 源匯項(xiàng)處理及水文地質(zhì)參數(shù)設(shè)置

模擬區(qū)地下水主要接受大氣降水和田間灌溉水補(bǔ)給，主要排泄項(xiàng)為蒸發(fā)和開(kāi)采。

(1)降雨入滲：廠址區(qū)年均降雨量為62.9 mm，降水入滲補(bǔ)給系數(shù)根據(jù)地表巖性參考水文地質(zhì)手冊(cè)推薦的經(jīng)驗(yàn)值。

(2)蒸發(fā)：由地下水位觀測(cè)資料，水位埋深大于極限蒸發(fā)深度4 m，故不考慮飽和地下水的蒸發(fā)。

(3)灌溉和開(kāi)采：模擬區(qū)內(nèi)180眼農(nóng)業(yè)灌溉井的地下水年開(kāi)采量約4 360萬(wàn)m3，計(jì)算得單井開(kāi)采量約664 m3/d，單井灌溉面積約199 800 m2，灌溉入滲系數(shù)參考水文地質(zhì)手冊(cè)的推薦值，取0.2，灌溉補(bǔ)給速率為(664 m3/d×0.2)÷199 800 m2=6.64×10-4m/d。即開(kāi)采量取664(m3/d)/眼，灌溉補(bǔ)給速率取6.64×10-4m/d。

(4)滲透系數(shù)：滲透系數(shù)根據(jù)水文地質(zhì)調(diào)查及水文地質(zhì)試驗(yàn)成果，結(jié)合含水層巖性，同時(shí)參考經(jīng)驗(yàn)值給出，并在模型識(shí)別驗(yàn)證中進(jìn)行微調(diào)。

4.3 模型識(shí)別和驗(yàn)證

建立地下水流模型后，通過(guò)擬合-校正進(jìn)行模型識(shí)別，本次模型識(shí)別與校正標(biāo)準(zhǔn)包括以下三個(gè)方面。

4.3.1識(shí)別后的模型概化和參數(shù)應(yīng)符合實(shí)際的水文地質(zhì)條件

校核出的含水層滲透系數(shù)和降雨入滲補(bǔ)給量分區(qū)如圖8和圖9所示。綜合考慮以往研究成果和現(xiàn)場(chǎng)試驗(yàn)數(shù)據(jù)及經(jīng)驗(yàn)值，各水文地質(zhì)參數(shù)基本符合實(shí)際特征。

圖8 滲透系數(shù)取值分區(qū)圖

圖9 降雨入滲補(bǔ)給量分區(qū)圖

4.3.2穩(wěn)定流狀態(tài)下觀測(cè)井地下水水位與計(jì)算水位基本一致

本次模擬擬合9眼長(zhǎng)期觀測(cè)孔，觀測(cè)地下水水位取一個(gè)水文年水位平均值，經(jīng)模型識(shí)別后的觀測(cè)點(diǎn)水位計(jì)算值和實(shí)測(cè)值擬合結(jié)果見(jiàn)圖10。圖中顯示實(shí)測(cè)水位與計(jì)算水位基本一致。

圖10 觀測(cè)孔水位擬合圖

4.3.3模擬地下水流場(chǎng)與調(diào)查得到的模擬區(qū)地下水流場(chǎng)基本一致

將模擬的地下水流場(chǎng)和實(shí)測(cè)的地下水等值線進(jìn)行對(duì)比，結(jié)果示于圖11。圖11中黑色線為根據(jù)實(shí)測(cè)水位數(shù)據(jù)繪制的地下水等值線，彩色線表示模擬水位分布。從對(duì)比圖可知，地下水模擬流場(chǎng)與實(shí)測(cè)流場(chǎng)形態(tài)基本一致，地下水流向相同。

圖11 地下水等水位線對(duì)比圖

5 包氣帶-飽和帶耦合模擬

5.1 核素遷移概念模型

核素遷移概念模型可概化為放射性液體通過(guò)包氣帶垂直進(jìn)入飽和地下水環(huán)境后，在對(duì)流、彌散、吸附和衰變作用下，隨地下水向環(huán)境敏感點(diǎn)遷移，核素遷移概念模型如圖12所示。

圖12 核素遷移概念模型示意圖

5.2 包氣帶一維垂向遷移模型與飽和帶三維模型的耦合

將包氣帶模擬獲得的潛水面處水流速度及C-14濃度隨時(shí)間變化值作為上邊界條件輸入到GMS飽和帶模型。

5.3 包氣帶-飽和帶耦合模擬結(jié)果及分析

C-14模擬結(jié)果示于圖13。由圖13可知C-14在對(duì)流作用下不斷向下游遷移，在彌散作用影響下，其污染范圍逐漸增大，在地下水稀釋及衰變、吸附等作用下，C-14活度濃度逐漸降低。事故發(fā)生后第1、5和10年，C-14放射性液體前緣遷移距離分別為1.5 m、4 m和4.2 m。

活度濃度單位為Bq/L，圖中網(wǎng)格非模型網(wǎng)格，間距1 m。

根據(jù)模擬數(shù)據(jù)結(jié)果，C-14在事故工況下發(fā)生泄漏后經(jīng)過(guò)包氣帶滲入飽和地下水環(huán)境，在地下水稀釋及吸附、衰變、對(duì)流、彌散等綜合作用下，其活度濃度大大降低，地下水中C-14的峰值活度濃度為2.15 Bq/L，對(duì)應(yīng)的年食入待積有效劑量為9.1×10-4mSv，遠(yuǎn)小于劑量限值0.1 mSv。事故工況后，放射性液體C-14活度濃度時(shí)空變化規(guī)律列于表3。

表3 事故工況后放射性液體C-14活度濃度時(shí)空變化規(guī)律

6 不考慮包氣帶的模擬研究

為了進(jìn)一步說(shuō)明包氣帶對(duì)核素遷移的遲滯作用，開(kāi)展不考慮包氣帶的模擬研究(以下簡(jiǎn)稱(chēng)飽水帶模擬研究)，假設(shè)放射性液體瞬時(shí)進(jìn)入飽和地下水環(huán)境，并將模擬結(jié)果與包氣帶-飽和帶耦合模擬結(jié)果進(jìn)行對(duì)比分析。

6.1 模擬條件

飽水帶模型與耦合模型的區(qū)別在于放射性液體進(jìn)入飽和地下水中的方式：耦合模型中放射性液體經(jīng)過(guò)包氣帶緩慢滲入飽和地下水，飽水帶模型中放射性液體瞬時(shí)進(jìn)入飽和地下水，亦即核素遷移模型的初始和源匯項(xiàng)(上邊界條件)設(shè)置。

6.1.1核素遷移模型初始條件

耦合模型中C-14初始活度濃度為0 Bq/L。而飽水帶模型中C-14初始活度濃度為9.4×104Bq/L，體積為8.50 m3。

6.1.2核素遷移模型源匯項(xiàng)

耦合模型中C-14補(bǔ)給活度濃度隨時(shí)間變化，詳見(jiàn)表3。而飽水帶模型將C-14活度濃度作為初始條件，不再考慮其補(bǔ)給源項(xiàng)。

6.2 模擬結(jié)果及分析

模擬結(jié)果示于圖14。由圖14和模擬數(shù)據(jù)結(jié)果可知：包氣帶對(duì)放射性核素的遷移存在明顯的遲滯作用，在包氣帶的影響下，地下水中C-14的活度濃度大大降低，放射性液體的遷移距離也明顯變小。

活度濃度單位為Bq/L，圖中網(wǎng)格非模型網(wǎng)格，間距10 m。

(1)在考慮包氣帶時(shí)，地下水中C-14峰值活度濃度為2.15 Bq/L，對(duì)應(yīng)的公眾食入年待積有效劑量為9.1×10-4mSv，小于評(píng)價(jià)準(zhǔn)則0.1 mSv；在不考慮包氣帶時(shí)，地下水中C-14峰值活度濃度為9.4×104Bq/L，對(duì)應(yīng)的年食入待積有效劑量為39.8 mSv，遠(yuǎn)大于0.1 mSv。

(2)在考慮包氣帶時(shí)，事故發(fā)生后第1、5和10年，C-14放射性液體遷移距離分別約1.5 m、4 m和4.2 m；在不考慮包氣帶時(shí)，事故發(fā)生后第1、5和10年，C-14放射性液體遷移距離分別約3 m、12 m和17 m。詳見(jiàn)表4。

表4 包氣帶對(duì)C-14活度濃度和遷移距離的影響

7 結(jié)論

為了能夠更加準(zhǔn)確地開(kāi)展地下水環(huán)境影響評(píng)價(jià)工作，本文采用擬三維數(shù)值方法開(kāi)展包氣帶-飽和帶遷移耦合模擬模型研究工作：包氣帶水分運(yùn)動(dòng)處理為一維垂向運(yùn)動(dòng)，飽和帶地下水為三維運(yùn)動(dòng)。為了進(jìn)一步說(shuō)明包氣帶對(duì)核素遷移的遲滯作用，開(kāi)展了不考慮包氣帶的模擬研究，并將模擬結(jié)果與包氣帶-飽和帶耦合模擬結(jié)果進(jìn)行對(duì)比分析。

主要結(jié)論如下：

(1)包氣帶模擬研究表明，由于包氣帶對(duì)核素遷移的遲滯作用(延遲和滯留)，對(duì)于吸附性較強(qiáng)的核素(Se-79、Pd-107、Cs-135和Sn-126)，其遷移至潛水面處的核素活度濃度遠(yuǎn)小于泄漏的放射性液體核素活度濃度，對(duì)應(yīng)的公眾食入年待積有效劑量已小于劑量限值0.1 mSv，只有吸附性較弱的C-14在潛水面處的核素活度濃度對(duì)應(yīng)年食入待積有效劑量大于0.1 mSv，需對(duì)其開(kāi)展耦合模擬研究工作。

(2)耦合模擬研究表明，C-14在對(duì)流作用下不斷向下游遷移，在彌散作用影響下，其污染范圍逐漸增大，在地下水稀釋及衰變、吸附等作用下，C-14活度濃度逐漸降低。地下水中C-14的峰值活度濃度為2.15 Bq/L，對(duì)應(yīng)的年食入待積有效劑量為9.1×10-4mSv，遠(yuǎn)小于劑量限值0.1 mSv。事故發(fā)生后第1、2和10年，C-14放射性液體前緣遷移距離分別為1.5 m、4 m和4.2 m。

(3)飽水帶模擬結(jié)果表明，包氣帶對(duì)放射性核素的遷移存在明顯的遲滯作用，在包氣帶的影響下，地下水中C-14的活度濃度大大降低，放射性液體的遷移距離也明顯變小：在考慮包氣帶時(shí)，地下水中C-14峰值活度濃度對(duì)應(yīng)的公眾食入年待積有效劑量為9.1×10-4mSv，小于評(píng)價(jià)準(zhǔn)則0.1 mSv，而不考慮包氣帶時(shí)，地下水中C-14峰值活度濃度對(duì)應(yīng)的年食入待積有效劑量為39.8 mSv，遠(yuǎn)大于0.1 mSv；在考慮包氣帶時(shí)，事故發(fā)生后第1、5和10年，C-14放射性液體遷移距離分別約1.5 m、4 m和4.2 m，而不考慮包氣帶時(shí)，事故發(fā)生后第1、5和10年，C-14放射性液體遷移距離分別約3 m、12 m和17 m。