凡倩瑩,劉 洋,馮永林,周繼強,王全榮*
(1.中國地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院,湖北 武漢 430078;2.山東省地質(zhì)科學(xué)研究院,山東 濟南 250013;3.自然資源部高寒干旱區(qū)礦山地質(zhì)環(huán)境修復(fù)工程技術(shù)創(chuàng)新中心,甘肅 蘭州 730000)
潛流帶是氮削減和河流水體凈化的主要場所,其水文地球化學(xué)過程是影響地下水環(huán)境演化過程的主控因素。近年來,工農(nóng)業(yè)活動產(chǎn)生大量的氮污染物[1],導(dǎo)致河流氮負(fù)荷逐年增加,河流水體富營養(yǎng)化[2-3]。河流下方及附近沉積物中河水與地下水的混合區(qū)(潛流帶)可以有效地降低或去除水環(huán)境中過量的氮[4-6]。自然水文或人為事件如洪水、降雨、潮汐作用、大壩釋水等引起的河水位瞬時改變[7],又進(jìn)一步增強了河水與地下水的混合,加劇了潛流帶溶質(zhì)的生物地球化學(xué)反應(yīng)。
高頻率的河水位波動會促進(jìn)河流與潛流帶間水量[8]、溶質(zhì)的交換[9]和遷移轉(zhuǎn)化[10-14],使得物質(zhì)濃度和反應(yīng)速率[15]在短距離和時間尺度內(nèi)發(fā)生變化,進(jìn)而影響硝酸鹽的去除[15-17]。數(shù)值模擬作為水文地質(zhì)學(xué)研究領(lǐng)域中常用的研究手段,因具有有效、快捷、靈活、相對廉價性等優(yōu)點,已被諸多學(xué)者用于研究河水位波動對潛流帶生物地球化學(xué)過程的影響。Gu等[16]首先采用數(shù)值模擬手段研究了河水位波動對潛流帶氮循環(huán)的影響,但研究中忽略了沉積物中原位碳釋放過程產(chǎn)生的溶解有機碳,也沒有考慮硝化作用;Shuai等[17]在Gu等研究的基礎(chǔ)上考慮了原位碳釋放過程和硝化作用,揭示了大壩釋水時河水位波動條件下潛流帶氮的遷移規(guī)律,但卻忽略了河床和含水層非均質(zhì)性在潛流帶氮循環(huán)中起到的影響,且研究中假設(shè)含水層是均質(zhì)的。當(dāng)前已有大量的野外觀測和研究表明,河床和含水層的非均質(zhì)性會影響潛流交換[18-19]、溶質(zhì)停留時間[20-22]和氮的去除程度[23-25]。一般來說,在非均質(zhì)場中滲透系數(shù)較大的區(qū)域,溶質(zhì)滯留時間較短且氧氣濃度較高[26],不利于反硝化作用的發(fā)生;滲透系數(shù)較小的區(qū)域,溶質(zhì)停留時間較長,有利于消耗氧氣形成厭氧環(huán)境,利于硝酸鹽的去除[27-29]。Sawyer等[30]采用數(shù)值模擬手段研究發(fā)現(xiàn),河床沉積物非均質(zhì)性影響進(jìn)入潛流帶內(nèi)的水流運動軌跡和在潛流帶內(nèi)的停留時間;Bardini等[31]采用數(shù)值模擬方法在研究河床形態(tài)對氮素遷移轉(zhuǎn)化的影響時發(fā)現(xiàn),泥沙滲透性和河流流速直接影響了養(yǎng)分的供應(yīng)和在河床中的停留時間,控制了沙丘下的反應(yīng)速率;Pescimoro等[32]在研究地梯度蜿蜒河流中具有砂和淤泥二元分布沉積物中的潛流交換時發(fā)現(xiàn),雖然更高淤泥占比的河床表現(xiàn)出較低的潛流交換率,但是沿流動路徑上硝酸鹽的去除效率更高;張強[33]建立了半切割河床式潛流帶模型,并將均質(zhì)潛流帶模型和低非均質(zhì)性潛流帶模型進(jìn)行了對比,發(fā)現(xiàn)潛流帶非均質(zhì)性降低了潛流交換量和硝酸鹽削減率;Wallace等[34]采用數(shù)值模擬手段研究了沉積物非均質(zhì)性對潮汐河流附近氮動態(tài)的影響,發(fā)現(xiàn)隨著含砂量的增加,河水中去除的硝酸鹽數(shù)量增加。
盡管目前關(guān)于河水位波動和滲透系數(shù)非均質(zhì)對潛流帶生物地球化學(xué)過程影響的研究已有很多,但大多研究往往忽略了河水位波動過程中氧氣溶解和原位碳釋放溶解有機碳的過程,這可能會導(dǎo)致非均質(zhì)潛流帶模型低估了潛流帶中溶質(zhì)供給和反硝化作用。為了進(jìn)一步探究潛流帶氮循環(huán)過程,本研究建立潛流帶二維流動和反應(yīng)傳輸模型,探究河水位波動下滲透系數(shù)非均質(zhì)性與潛流帶內(nèi)氮素遷移轉(zhuǎn)化過程的空間關(guān)系,了解滲透系數(shù)非均質(zhì)性對潛流交換的影響、溶質(zhì)在空間上遷移轉(zhuǎn)化的差異和硝酸鹽去除率的變化。
使用COMSOL Multiphysics構(gòu)建二維潛流帶概念模型,用于模擬河流波動引起的潛流帶內(nèi)氮素遷移轉(zhuǎn)化過程。模型幾何形狀取自Hornsby Bend站點上用全站儀測量的從河漫灘到河道中心的河岸高程輪廓,具體見圖1。
注:透水層外部高度為河水位高度;開放邊界外部濃度為河水濃度和以氣泡形式捕獲的空氣濃度。圖1 二維潛流帶氮循環(huán)過程示意圖(改自Shuai等[17])Fig.1 Schematic diagram of the nitrogen cycle process in a two-dimensional hyporheic zone (adapted from Shuai et al.[17])
本研究采用Richard方程和對流-彌散-反應(yīng)方程刻畫潛流帶內(nèi)水流和溶質(zhì)的遷移過程,其控制方程如下:
(1)
(2)
式中:p為壓力(Pa);g為重力加速度(m/s2);ρ為流體密度(kg/m3);μ為流體動力黏度(Pa·s);Cm為比水容量(1/m);Se為有效飽和度(無量綱);S為儲水系數(shù)(1/Pa);ks為飽和水力滲透率(m2);kr為相對滲透率(無量綱),是有效飽和度的函數(shù);t為時間(s);z為位置高度(m);Qm為源匯項[kg/(m3·s)];θ為含水率(m3/m3);Ci為物質(zhì)i的濃度(mg/L);u為達(dá)西流速(m/s);D為水動力彌散系數(shù)(m2/s);Ri為物質(zhì)i的反應(yīng)速率[mg/(L·d)]。
有氧呼吸作用 CH2O+O2→CO2+H2O
原位有機物或總有機碳的溶解可能導(dǎo)致DOC的輸入,該過程通過一級傳質(zhì)過程來模擬[17,37-40]:
(3)
式中:CDOC和CTOC分別為沉積物中DOC和總有機碳(TOC)的濃度(mg/L);α為一階傳質(zhì)系數(shù)(1/h);Kd為DOC分布系數(shù)(L/kg)。
有氧呼吸作用、硝化作用和反硝化作用由多重Monod動力學(xué)描述[41]。存在DO時,使用非競爭性來抑制反硝化作用[42]。此外,系統(tǒng)中的氧通常用于能量優(yōu)先級更高的反應(yīng),如有氧呼吸;而剩余的氧則作為具有較低能量優(yōu)先級的二級溶質(zhì)被消耗,如反硝化反應(yīng)。而溶解氧需求的分布系數(shù)yO2是基于氧化和硝化反應(yīng)的吉布斯自由能計算的。綜上所述,各化學(xué)反應(yīng)項表示如下:
(4)
(5)
(6)
(7)
基于概念模型和數(shù)學(xué)模型,首先將潛流帶模型離散化。為了提高模型模擬精度和減少計算量,采用三角形網(wǎng)格進(jìn)行剖分(圖2),上邊界附近最小空間步長為0.1 m,遠(yuǎn)離上邊界最大空間步長為0.3 m,將模型離散為10 184個網(wǎng)格。模型模擬總時長為24 h,時間步長為1 h。
對于水流模型,河道中心剖面與底部邊界處理為隔水邊界;右側(cè)邊界和上邊界分別設(shè)置為定水頭邊界和透水層邊界,不考慮降雨或蒸發(fā)情況(具體設(shè)置參考Cardenas等[43])。河水位波動近似表示成振幅為0.5 m、周期為24 h的正弦函數(shù)。
對于溶質(zhì)運移模型,河道中心剖面與底部邊界處理為零通量邊界;右側(cè)邊界設(shè)置為流出邊界;上邊界設(shè)置為開放邊界,在這里允許發(fā)生溶質(zhì)的流入和流出,且外部濃度等于河流中溶質(zhì)濃度。假定在地下水水位上升時,與大氣相連的非飽和帶中的DO處于飽和狀態(tài)。
表1 水流場和溶質(zhì)運移場主要參數(shù)設(shè)置
滲透系數(shù)的非均質(zhì)性影響河流-潛流帶的水動力交換和流場分布,一個24 h周期內(nèi)河水位波動和潛流交換通量,見圖3。
注:正值表示含水層向河流補給(出滲),負(fù)值表示河水向含水層補給(入滲);平均河水位和零通量則分別用紅色和藍(lán)色虛線表示。圖3 一個24 h周期內(nèi)河水位波動(藍(lán)色實線)和潛流 交換通量(點線)Fig.3 Fluctuations in river level (solid blue line) and hyporheic exchange flux (dotted line) in a 24-hour period
由圖3可知:河水位上升時,河水向含水層入滲(負(fù)值),河水位下降至平均河水位以下后,水由河水入滲含水層(負(fù)值)轉(zhuǎn)變?yōu)榱鞒龊畬?正值)向河水補給,并在18 h時河水位達(dá)到最低;18 h后,河水位雖然有所上升,但此時河水位仍低于平均河水位,含水層入滲補給河水;河流-含水層交界面上的入滲和出滲通量隨著河水位波動而上升和下降,在河段河水位達(dá)到最大值前入滲流量達(dá)到峰值,而出滲通量在河段河水位達(dá)到最小值前達(dá)到峰值,這一趨勢與 Musial等[9]受潮汐影響的近溪流實地觀察結(jié)果和Shuai等[17]的模擬結(jié)果一致。
注:黑色箭頭為流線。圖4 不同方差下潛流帶滲透系數(shù)場和流場分布圖Fig.4 Diagram of permeability coefficient fields and flow fields of hyporheic zones under different variance values
由圖4可知,隨著潛流帶含水層非均質(zhì)性的增強,流線越曲折,此外,含水層的非均質(zhì)性還影響潛流帶的交換速率,含水層的非均質(zhì)性越強,潛流帶交界面處的等水位線向兩側(cè)延申越遠(yuǎn),潛流帶交換范圍越廣,反應(yīng)越劇烈。
圖5 不同方差下潛流帶中溶質(zhì)運移過程圖Fig.5 Solute transport process in hyporheic zones with different variance values
圖6 不同方差下潛流帶中觀測點處溶質(zhì)濃度和化學(xué)反應(yīng)速率對比圖Fig.6 Comparison between solute concentration and chemical reaction rate at observation points in hyporheic zone under different variance values
硝化作用受溶解氧控制,主要發(fā)生在河流-含水層界面附近的富氧區(qū)域[圖5(d)];反硝化作用發(fā)生在厭氧區(qū)域,主要沿著有氧呼吸區(qū)域外圍以帶狀的鋒面形式分布[圖5(e)],越遠(yuǎn)離高氧區(qū),反硝化反應(yīng)速率就越高。滲透系數(shù)的非均質(zhì)性使?jié)摿鲙?nèi)高滲透層增多,溶質(zhì)的滲透深度增加,氨氮、硝酸鹽等物質(zhì)向潛流帶深處運移,這在一定程度上了改變了硝化和反硝化帶的形態(tài),擴大了反應(yīng)的作用區(qū)。與潛流帶內(nèi)溶質(zhì)濃度的空間分布形狀相類似,隨著含水層非均質(zhì)性的增強,潛流帶內(nèi)高滲透層增加,硝化作用和反硝化作用區(qū)分布形狀越不規(guī)則,邊緣越“曲折”。
為了定量分析潛流帶對硝酸鹽的去除情況,Shuai等[17]假定硝酸鹽去除效率是含水層中去除的硝酸鹽總量與進(jìn)入含水層中硝酸鹽總量的比值NRE,其可表示為
其中,反硝化作用去除的含水層中硝酸鹽量MDN和經(jīng)河流-含水層界面進(jìn)入含水層中的硝酸鹽總量Min分別表示如下:
滲透系數(shù)的非均質(zhì)性影響了潛流帶內(nèi)硝酸鹽的去除,不同方差下潛流帶內(nèi)硝酸鹽的去除情況見表2。
表2 不同方差下潛流帶內(nèi)硝酸鹽去除情況
由表2可知:隨著含水層非均質(zhì)性程度的增加,硝化、反硝化反應(yīng)速率和溶質(zhì)滲透深度均有所增加??偟脕碚f,一個周期內(nèi)河流向潛流帶輸入的硝酸鹽總量和經(jīng)反硝化作用去除的硝酸鹽總量也隨之增加,潛流帶內(nèi)硝酸鹽去除率隨著滲透系數(shù)非均質(zhì)性程度的增加而增大。
本研究利用數(shù)值模擬手段研究了滲透系數(shù)非均質(zhì)性對潛流帶氮素遷移轉(zhuǎn)化和硝酸鹽去除的影響。主要得到結(jié)論如下:
1) 潛流內(nèi)滲透系數(shù)非均質(zhì)性程度的增加加快了潛流交換并提供了更高的反應(yīng)速率。隨著非均質(zhì)性程度的增加,潛流帶滲流面上入滲率和出滲率也越大,河水-含水層間的水力交換越劇烈。
2) 潛流內(nèi)滲透系數(shù)非均質(zhì)性影響了溶質(zhì)在含水層中的運移情況,控制了潛流帶內(nèi)的反應(yīng)動力學(xué),改變了溶質(zhì)分布形狀和硝化、反硝化作用的發(fā)生區(qū)域,隨著非均質(zhì)性的加強,溶質(zhì)濃度和反應(yīng)區(qū)在空間的分布越不規(guī)則,分布區(qū)邊緣越“曲折”。
3) 滲透系數(shù)的非均質(zhì)性影響潛流帶去除硝酸鹽的能力,更高的滲透系數(shù)非均質(zhì)性表現(xiàn)出更好的硝酸鹽去除效果,硝酸鹽去除率也有所增加。