放射性核素遷移的包氣帶-飽和帶耦合模擬研究

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

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

為保證核設(shè)施建設(shè)、運(yùn)行及退役等階段土壤及地下水環(huán)境的安全，根據(jù)國家核安全法規(guī)、標(biāo)準(zhǔn)和導(dǎo)則的要求，必須通過實(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ì)于放射性核素來說，包氣帶對(duì)放射性核素的阻滯作用非常明顯，李書紳等[1]的研究表明，在自然條件下，可吸附核素的遷移距離非常短，因此研究放射性核素在包氣帶中遷移是必要的。放射性核素在包氣帶中的遷移主要以垂直向?yàn)橹?，且與地下水環(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)的問題進(jìn)行耦合而得到擬三維數(shù)值模型。包氣帶一維垂向模擬選用Hydrus-1D軟件，飽和帶三維模擬選用GMS軟件。將潛水面作為包氣帶模型下邊界，同時(shí)作為飽和帶模型上邊界，通過包氣帶模擬獲得包氣帶下邊界地下水流量及放射性核素濃度隨時(shí)間變化值，并將其作為飽和帶的源匯項(xiàng)(上邊界條件)輸入到飽和帶模型，實(shí)現(xiàn)包氣帶-飽和帶的松散耦合。

1 國內(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]采用類似的方法建立了一個(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)簡化為一維垂向運(yùn)動(dòng)，保持飽和帶三維運(yùn)動(dòng)的特點(diǎn)，將兩者進(jìn)行耦合。查元源(2014)[11]將Richards方程模型及水均衡模型與三維地下水模型耦合，建立了區(qū)域飽和-非飽和三維水分運(yùn)動(dòng)耦合模型。于洋等(2017)[12]通過GMS和Hydrus的結(jié)合，實(shí)現(xiàn)了包氣帶-飽和帶水鹽運(yùn)移的耦合模擬。張雅楠(2019)[13]明晰了HYDRUS-MODFLOW模型的耦合原理及運(yùn)行機(jī)制。

2 研究區(qū)概況

研究區(qū)屬大陸性干旱氣候，降水少，蒸發(fā)強(qiáng)烈。地勢北高南低，北部為基巖山地，海拔1 250.0～1 440.0 m，南部為碎石戈壁，海拔一般為1 229.0～1 275.0 m。多年平均降雨量62.9 mm，多年平均蒸發(fā)量為2 338.9 mm。地下水可分為松散巖類孔隙水和基巖裂隙水兩類。松散巖類孔隙水區(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)場地包氣帶厚度約3 m，假設(shè)罐體中的放射性液體發(fā)生泄漏，并由地表垂直入滲補(bǔ)給包氣帶，泄漏總量為8.50 m3，核素種類包括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通過現(xiàn)場土-水特征曲線試驗(yàn)和現(xiàn)場包氣帶滲水試驗(yàn)獲取初值，并利用包氣帶水分遷移試驗(yàn)，通過模型識(shí)別驗(yàn)證獲取模型參數(shù)取值；干密度ρd取測試值；彌散度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é)果，接下來只對(duì)潛水面處核素活度濃度對(duì)應(yīng)年待積有效劑量大于評(píng)價(jià)準(zhǔn)則(0.1 mSv)的C-14開展耦合模擬研究。

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)沒有明顯的隔水層，根據(jù)地下水位長期監(jiān)測結(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ā)和開采。

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

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

(3)灌溉和開采：模擬區(qū)內(nèi)180眼農(nóng)業(yè)灌溉井的地下水年開采量約4 360萬m3，計(jì)算得單井開采量約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。即開采量取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)證

建立地下水流模型后，通過擬合-校正進(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)場試驗(yàn)數(shù)據(jù)及經(jīng)驗(yàn)值，各水文地質(zhì)參數(shù)基本符合實(shí)際特征。

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

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

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

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

圖10 觀測孔水位擬合圖

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

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

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

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

5.1 核素遷移概念模型

核素遷移概念模型可概化為放射性液體通過包氣帶垂直進(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)過包氣帶滲入飽和地下水環(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)一步說明包氣帶對(duì)核素遷移的遲滯作用，開展不考慮包氣帶的模擬研究(以下簡稱飽水帶模擬研究)，假設(shè)放射性液體瞬時(shí)進(jìn)入飽和地下水環(huán)境，并將模擬結(jié)果與包氣帶-飽和帶耦合模擬結(jié)果進(jìn)行對(duì)比分析。

6.1 模擬條件

飽水帶模型與耦合模型的區(qū)別在于放射性液體進(jìn)入飽和地下水中的方式：耦合模型中放射性液體經(jīng)過包氣帶緩慢滲入飽和地下水，飽水帶模型中放射性液體瞬時(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í)間變化，詳見表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。詳見表4。

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

7 結(jié)論

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

(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。