氡示蹤河底地下水-河水相互交換數值模擬研究

＊樂慶國

（香港城市大學 香港 999077）

1.概念模型

(1)河水-地下水交換模型

天然河流里，具有周期的三角沙丘河床地形相當普遍。以河流流向為X軸，垂直河床面為Y軸，建立二維河床剖面模型。模型中，擬定河床沉積物中存在連續的三角形沙丘微地形，緩坡方向與河流流向相同，陡坡方向與河流流向相反。邊界條件中，以模型左邊界為速度入口，模型右邊界為壓力出口，河床、河面為無滑移邊界，河床由5個連續的三角形沙丘構成，每個沙丘波長1m，位于0.75m處有一尖峰，高為0.115m。采用ANSYS Fluent模擬河流紊流流動，采用的剖分方法為三角剖分，節點數在90000個左右，盡可能避免由于網格剖分數量給模擬結果帶來影響，以此模型模擬計算河水給其下伏的潛流帶的壓力及其分布。

(2)河流紊流模型

ANSYS Fluent模擬河流紊流原理為有限體積法求解不可壓縮流體納維-斯托克斯方程（Navier-Stocks Equation）。所用到的方程為RANS方程和k-ω湍流模型，RANS方程如下：

式中，ρ、μ為流體的密度和動力黏度；Ui（i=1,2）、ui（i=1,2）分別為流體的平均速度和瞬時速度；g為重力加速度；P為流體平均壓強；Sij為應變速率張量，其式為：

RANS方程與k湍流動能和湍流耗散率ω有關，其式為：

其中，τij為雷諾切應力；vt為渦流運動黏度；δij為克羅內克符號。

湍流動能κ方程為：

湍流耗散率ω方程為：

其中，μt=ρvt，α=5/9，β=0.075，β*=0.09，σκ=σω=0.5。

(3)水流和溶質運移模型

河床形態引起的潛流交換主要包括泵吸交換（Pumping Exchange）和沖淤交換（Turnover Exchange）。假設三角沙丘河床不發生移動，潛流交換主要由泵吸作用引起。孔隙介質飽和滲流可用方程描述：

假設河床沉積物為均質各向同性，h為水頭；K為滲透系數；n為孔隙度。

222Rn在孔隙介質中運移涉及對流和彌散，且由于其具有放射性，它的產生和衰變方程為：

2.河流流速對氡示蹤潛流交換的影響

(1)模擬方案

水動力條件作為潛流帶內物質和能量交換的重要驅動力，不僅控制著地表水和地下水的水量交換，影響兩者的混合比例，同時控制著隨水流進入潛流帶內各物質的停留時間、混合運移和時空分布模式。在本研究中，河床面初始壓力與河流流速直接相關聯，不同的河流流速會產生不同的河床面初始壓力，初始壓力在COMSOL Multiphysics中以壓力邊界的形式存在，壓力邊界是邊界條件的一種，是確定一個水文地質模型能否求解的一個關鍵參數。本研究選定河流流速為變量，定量研究河流流速對初始壓力的影響，并探究后者對于222Rn示蹤潛流交換的影響。模擬方案的取值見表1，其余參數恒定不變，見表2。

表1 各個模型河流流速取值

表2 其余參數取值

(2)結果與分析

圖1 不同河流流速下河床的壓力分布

①河流流速變化下河床壓力分布

從圖2來看，四個模型都展現出了河水在一定流速驅動下，由于沙丘形態，流速橫截面縮小，沖擊著迎水面的沙丘，在大約沙丘迎水面的中部壓力最大，河水由此進入河床沉積物，在沙丘的左端和右端壓力最小，河水在此返回河流。沙丘迎水面和背水面的壓力差隨著流速的增加而呈平方增長每一個沙丘的壓力線型十分相似，但是在總體趨勢上呈現出壓力隨著河水流動距離增加而減少，河水每流經一個沙丘，沙丘壓力峰值下降大約20%，最小值則下降大約50%。值得注意的是，在第三個沙丘的右側0.75m處低凹區，四個模型都出現了負壓。此外，模型1，2在沙丘頂0.75m處都出現了顯著低于左右側的壓力值，模型3，4則不明顯或基本可以忽略。河床面的壓力值與河流流速的平方呈正相關關系，模型2和模型4的河水流速只相差一個量級，但河床壓力數值卻相差兩個量級，這意味著河流流速增加導致的河床面壓力增加的邊際效應十分顯著。

圖2 不同河流流速下潛流帶的水頭分布圖

由于不同河流流速下每個沙丘的壓力線型相對一致，四個模型的水頭分布圖較為相似，具體表現為：每個沙丘中部0.5m處水頭較大，作為系統的源；沙丘左右側0～0.1m和0.9～1m水頭較小，作為系統的匯；較淺的區域水頭差較大，隨著沉積物深度增加，水力梯度逐漸減小，潛流交換強度隨著沉積物深度增加而降低。由于模型3，4的第三個沙丘的低壓力區相較于前兩個沙丘壓力差過于巨大，低壓力區已蔓延到模型底部，意味著來自前兩個沙丘的河水更有可能會深入潛流帶底部與深部地下水進行交換，再由模型右側的低壓力區離開。水頭數值與初始壓力呈正相關關系。

②河流流速的變化對氡示蹤潛流交換的影響

圖4至圖6可以看出，其相同點是，在模擬的早期，222Rn示蹤效果較好，往往在模擬時間t=1d時，河水就已擴散到潛流帶中相當大的一部分，這部分區域可以用低濃度222Rn示蹤出來，即該區域的222Rn濃度遠低于地下水中的222Rn濃度，約1至2個量級；越往深處，由于含水層固體顆粒礦物（226Ra，222Rn的母體）衰變產生222Rn進入潛流帶，潛流帶的222Rn濃度逐漸升高，逐漸接近地下水氡濃度，氡示蹤潛流的交換的效果逐漸減弱。

模擬期第10天已基本展現出穩定的態勢，相較于第20天模擬期結束的最后結果相差不大，甚至完全一致；由上到下，從低濃度區到高濃度區，每一個模型都出現了一層漸變過渡帶，近似呈“3V”型，高濃度區的氡經由這一漸變帶擴散至低濃度區，并在此過程中濃度逐步降低。原因為222Rn經歷了數個半衰期，從模擬一開始時進入模型中的222Rn在抵達濃度漸變帶區域基本衰變殆盡，對于濃度漸變帶及其之下的深度，222Rn對于示蹤潛流交換的效果不佳。

針對不同區域，從圖3到圖6，可以清楚的發現初始壓力和氡濃度分布之間的關系：初始壓力越大，低濃度的222Rn區域就越大。模型1對應極低流速下的河流，例如，河流下游接近入海口的河段。由于水動力條件較為薄弱，在滲透系數不太大的情況下，河水向潛流帶擴散有限，僅能影響河床以下約0.1m區域，濃度漸變帶約0.3m，對于222Rn來說，示蹤再深處的區域則力不能及。模型2展現的氡濃度分布圖相較模型1，低222Rn濃度區域明顯大得多，大約0.2m寬，漸變帶寬度約0.4m。模型3，4結果則不相同，在模擬期結束的第20天，低222Rn濃度區域都擴散到了近模型底部，只是模型4的擴散速度更快，它們的低濃度區域和濃度漸變帶也較模型1，2更厚，約0.5m。

圖3 模型1河流流速為0.06m/s 222Rn濃度分布圖

圖4 模型2河流流速為0.2m/s 222Rn濃度分布圖

圖5 模型3河流流速為0.5m/s 222Rn濃度分布圖

圖6 模型4河流流速為1m/s 222Rn濃度分布圖

造成模擬結果如此分異的原因為河床面上的壓力差不同，模型1，2的沙丘迎水面—背水面壓力差較小，難以驅使河水深入潛流帶進行交換；而模型3，4的壓力差高達數百帕，這有助于提高潛流交換強度，以對流為主的溶質遷移速率也會加強，在濃度圖上的低濃度區域也相對較大。而在這種情況下，222Rn就能夠很好地示蹤模型中的潛流交換。

3.結論

基于COMSOL Multiphysics建立的河水與地下水水流-溶質運移耦合模型模擬結果，表明在模擬的早期，潛流帶立即形成一片由河水構成的低濃度區，隨著河水進一步向深處擴散，氡的運移減慢，低濃度區擴張逐漸放緩穩定，并從上至下形成一條濃度由低變高的漸變帶。ANSYS Fluent模擬河流紊流的結果表明，地表水流條件和微地形形態控制著初始壓力，其中初始壓力對于前者尤其敏感，與前者的平方呈正相關；后者則會改變初始壓力的分布。二者皆對氡示蹤潛流交換具有顯著影響，河水水流增大，河床地形坡度增大，氡示蹤效果相對更加顯著。

水文地質參數對于氡示蹤潛流交換的效果具有重要影響。滲透系數K≤0.1，介質滲透性較差，模型中低濃度區域僅局限于水沙界面，222Rn的示蹤效果不佳；如果K≥1，低濃度區域顯著大于前者，且滲透系數越大，低濃度區域越大，222Rn的示蹤效果越為明顯。孔隙度影響著多孔介質的放射性元素衰變產生的氡進入地下水的速率，孔隙度越大則氡的產生速率越大，越大的孔隙度會使氡的濃度平衡深度變淺。彌散系數控制著在進入含水層的河水在對流路徑上的彌散，彌散系數越大，河水中的氡能遷移的更遠、更深，低濃度區域也越大，222Rn的示蹤效果也會越好。地下水氡濃度值對于氡的示蹤影響甚微，其值越小，越難以與地表水區別開，222Rn的示蹤效果也越差。通過對以上地表水條件和水文地質參數變化影響的研究，刻畫了每個參數對于氡示蹤地表水—地下水的交換的影響，可為地下水與地表水的交叉研究和實踐提供理論基礎。