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

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

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

為保證核設施建設、運行及退役等階段土壤及地下水環境的安全，根據國家核安全法規、標準和導則的要求，必須通過實驗和數值模擬等手段，充分掌握核設施廠址地區放射性核素在土壤及水環境中的遷移擴散規律，提供工程設計所需的參數資料。

以往對假想事故條件下放射性液體的地下水環境影響評價，假設放射性液體瞬時進入飽和地下水環境，不考慮包氣帶對核素遷移的遲滯作用。對于放射性核素來說，包氣帶對放射性核素的阻滯作用非常明顯，李書紳等[1]的研究表明，在自然條件下，可吸附核素的遷移距離非常短，因此研究放射性核素在包氣帶中遷移是必要的。放射性核素在包氣帶中的遷移主要以垂直向為主，且與地下水環評的評價范圍相比，水平向的運動幾乎可以忽略。

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

1 國內外研究進展

大多數飽和-非飽和水流耦合模型都基于MODFLOW建立，基于一維非飽和水流運動和三維地下水運動的含水層數值模擬系統如圖1所示。

圖1 基于MODFLOW建立的擬三維飽和-非飽和水流系統

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

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

2 研究區概況

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

3 包氣帶模擬研究

3.1 數值模型的建立

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

試驗場地包氣帶厚度約3 m，假設罐體中的放射性液體發生泄漏，并由地表垂直入滲補給包氣帶，泄漏總量為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)假設初始條件下水流處于穩定狀態，水流模型初始值由壓力水頭給出：h(z,0)=h0(z)，-H≤z≤0。式中：h為土壤水壓力冰頭，H為包氣帶厚度，z為垂直坐標，設潛水面處壓力水頭為0，向上遞減，地表處壓力水頭為-3 m。

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

3.3 邊界條件

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

(1)

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

下邊界選擇潛水面，設為定水頭邊界：

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

(2)

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

(3)

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

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

3.4 模型主要參數確定

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

表1 參數取值表

表2 分配系數和衰變常數取值

3.5 模擬結果及分析

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

圖2 潛水面處水流速度隨時間變化

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

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

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

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

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

由圖2～7模擬數據結果可知：

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

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

即包氣帶對核素的遷移存在明顯的遲滯作用，根據包氣帶模擬結果，接下來只對潛水面處核素活度濃度對應年待積有效劑量大于評價準則(0.1 mSv)的C-14開展耦合模擬研究。

4 飽和帶地下水水流數值模型

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

模擬范圍與邊界條件充分考慮評價區水文地質條件確定。

模擬區內地下水運動符合達西定律。地下水系統的垂向運動主要是大氣降水在地表入滲后向下補給含水層，孔隙水含水層垂向補給基巖裂隙水，模擬區內沒有明顯的隔水層，根據地下水位長期監測結果，該地區地下水水位隨季節變化不明顯。因此，模擬區地下水系統可概化為非均質三維穩態地下水流系統。

4.2 源匯項處理及水文地質參數設置

模擬區地下水主要接受大氣降水和田間灌溉水補給，主要排泄項為蒸發和開采。

(1)降雨入滲：廠址區年均降雨量為62.9 mm，降水入滲補給系數根據地表巖性參考水文地質手冊推薦的經驗值。

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

(3)灌溉和開采：模擬區內180眼農業灌溉井的地下水年開采量約4 360萬m3，計算得單井開采量約664 m3/d，單井灌溉面積約199 800 m2，灌溉入滲系數參考水文地質手冊的推薦值，取0.2，灌溉補給速率為(664 m3/d×0.2)÷199 800 m2=6.64×10-4m/d。即開采量取664(m3/d)/眼，灌溉補給速率取6.64×10-4m/d。

(4)滲透系數：滲透系數根據水文地質調查及水文地質試驗成果，結合含水層巖性，同時參考經驗值給出，并在模型識別驗證中進行微調。

4.3 模型識別和驗證

建立地下水流模型后，通過擬合-校正進行模型識別，本次模型識別與校正標準包括以下三個方面。

4.3.1識別后的模型概化和參數應符合實際的水文地質條件

校核出的含水層滲透系數和降雨入滲補給量分區如圖8和圖9所示。綜合考慮以往研究成果和現場試驗數據及經驗值，各水文地質參數基本符合實際特征。

圖8 滲透系數取值分區圖

圖9 降雨入滲補給量分區圖

4.3.2穩定流狀態下觀測井地下水水位與計算水位基本一致

本次模擬擬合9眼長期觀測孔，觀測地下水水位取一個水文年水位平均值，經模型識別后的觀測點水位計算值和實測值擬合結果見圖10。圖中顯示實測水位與計算水位基本一致。

圖10 觀測孔水位擬合圖

4.3.3模擬地下水流場與調查得到的模擬區地下水流場基本一致

將模擬的地下水流場和實測的地下水等值線進行對比，結果示于圖11。圖11中黑色線為根據實測水位數據繪制的地下水等值線，彩色線表示模擬水位分布。從對比圖可知，地下水模擬流場與實測流場形態基本一致，地下水流向相同。

圖11 地下水等水位線對比圖

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

5.1 核素遷移概念模型

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

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

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

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

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

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

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

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

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

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

為了進一步說明包氣帶對核素遷移的遲滯作用，開展不考慮包氣帶的模擬研究(以下簡稱飽水帶模擬研究)，假設放射性液體瞬時進入飽和地下水環境，并將模擬結果與包氣帶-飽和帶耦合模擬結果進行對比分析。

6.1 模擬條件

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

6.1.1核素遷移模型初始條件

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

6.1.2核素遷移模型源匯項

耦合模型中C-14補給活度濃度隨時間變化，詳見表3。而飽水帶模型將C-14活度濃度作為初始條件，不再考慮其補給源項。

6.2 模擬結果及分析

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

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

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

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

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

7 結論

為了能夠更加準確地開展地下水環境影響評價工作，本文采用擬三維數值方法開展包氣帶-飽和帶遷移耦合模擬模型研究工作：包氣帶水分運動處理為一維垂向運動，飽和帶地下水為三維運動。為了進一步說明包氣帶對核素遷移的遲滯作用，開展了不考慮包氣帶的模擬研究，并將模擬結果與包氣帶-飽和帶耦合模擬結果進行對比分析。

主要結論如下：

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

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

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