凡倩瑩,劉 洋,馮永林,周繼強(qiáng),王全榮*

(1.中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院,湖北 武漢 430078;2.山東省地質(zhì)科學(xué)研究院,山東 濟(jì)南 250013;3.自然資源部高寒干旱區(qū)礦山地質(zhì)環(huán)境修復(fù)工程技術(shù)創(chuàng)新中心,甘肅 蘭州 730000)

潛流帶是氮削減和河流水體凈化的主要場(chǎng)所,其水文地球化學(xué)過程是影響地下水環(huán)境演化過程的主控因素。近年來,工農(nóng)業(yè)活動(dòng)產(chǎn)生大量的氮污染物[1],導(dǎo)致河流氮負(fù)荷逐年增加,河流水體富營(yíng)養(yǎng)化[2-3]。河流下方及附近沉積物中河水與地下水的混合區(qū)(潛流帶)可以有效地降低或去除水環(huán)境中過量的氮[4-6]。自然水文或人為事件如洪水、降雨、潮汐作用、大壩釋水等引起的河水位瞬時(shí)改變[7],又進(jìn)一步增強(qiáng)了河水與地下水的混合,加劇了潛流帶溶質(zhì)的生物地球化學(xué)反應(yīng)。

高頻率的河水位波動(dòng)會(huì)促進(jìn)河流與潛流帶間水量[8]、溶質(zhì)的交換[9]和遷移轉(zhuǎn)化[10-14],使得物質(zhì)濃度和反應(yīng)速率[15]在短距離和時(shí)間尺度內(nèi)發(fā)生變化,進(jìn)而影響硝酸鹽的去除[15-17]。數(shù)值模擬作為水文地質(zhì)學(xué)研究領(lǐng)域中常用的研究手段,因具有有效、快捷、靈活、相對(duì)廉價(jià)性等優(yōu)點(diǎn),已被諸多學(xué)者用于研究河水位波動(dòng)對(duì)潛流帶生物地球化學(xué)過程的影響。Gu等[16]首先采用數(shù)值模擬手段研究了河水位波動(dòng)對(duì)潛流帶氮循環(huán)的影響,但研究中忽略了沉積物中原位碳釋放過程產(chǎn)生的溶解有機(jī)碳,也沒有考慮硝化作用;Shuai等[17]在Gu等研究的基礎(chǔ)上考慮了原位碳釋放過程和硝化作用,揭示了大壩釋水時(shí)河水位波動(dòng)條件下潛流帶氮的遷移規(guī)律,但卻忽略了河床和含水層非均質(zhì)性在潛流帶氮循環(huán)中起到的影響,且研究中假設(shè)含水層是均質(zhì)的。當(dāng)前已有大量的野外觀測(cè)和研究表明,河床和含水層的非均質(zhì)性會(huì)影響潛流交換[18-19]、溶質(zhì)停留時(shí)間[20-22]和氮的去除程度[23-25]。一般來說,在非均質(zhì)場(chǎng)中滲透系數(shù)較大的區(qū)域,溶質(zhì)滯留時(shí)間較短且氧氣濃度較高[26],不利于反硝化作用的發(fā)生;滲透系數(shù)較小的區(qū)域,溶質(zhì)停留時(shí)間較長(zhǎng),有利于消耗氧氣形成厭氧環(huán)境,利于硝酸鹽的去除[27-29]。Sawyer等[30]采用數(shù)值模擬手段研究發(fā)現(xiàn),河床沉積物非均質(zhì)性影響進(jìn)入潛流帶內(nèi)的水流運(yùn)動(dòng)軌跡和在潛流帶內(nèi)的停留時(shí)間;Bardini等[31]采用數(shù)值模擬方法在研究河床形態(tài)對(duì)氮素遷移轉(zhuǎn)化的影響時(shí)發(fā)現(xiàn),泥沙滲透性和河流流速直接影響了養(yǎng)分的供應(yīng)和在河床中的停留時(shí)間,控制了沙丘下的反應(yīng)速率;Pescimoro等[32]在研究地梯度蜿蜒河流中具有砂和淤泥二元分布沉積物中的潛流交換時(shí)發(fā)現(xiàn),雖然更高淤泥占比的河床表現(xiàn)出較低的潛流交換率,但是沿流動(dòng)路徑上硝酸鹽的去除效率更高;張強(qiáng)[33]建立了半切割河床式潛流帶模型,并將均質(zhì)潛流帶模型和低非均質(zhì)性潛流帶模型進(jìn)行了對(duì)比,發(fā)現(xiàn)潛流帶非均質(zhì)性降低了潛流交換量和硝酸鹽削減率;Wallace等[34]采用數(shù)值模擬手段研究了沉積物非均質(zhì)性對(duì)潮汐河流附近氮?jiǎng)討B(tài)的影響,發(fā)現(xiàn)隨著含砂量的增加,河水中去除的硝酸鹽數(shù)量增加。

盡管目前關(guān)于河水位波動(dòng)和滲透系數(shù)非均質(zhì)對(duì)潛流帶生物地球化學(xué)過程影響的研究已有很多,但大多研究往往忽略了河水位波動(dòng)過程中氧氣溶解和原位碳釋放溶解有機(jī)碳的過程,這可能會(huì)導(dǎo)致非均質(zhì)潛流帶模型低估了潛流帶中溶質(zhì)供給和反硝化作用。為了進(jìn)一步探究潛流帶氮循環(huán)過程,本研究建立潛流帶二維流動(dòng)和反應(yīng)傳輸模型,探究河水位波動(dòng)下滲透系數(shù)非均質(zhì)性與潛流帶內(nèi)氮素遷移轉(zhuǎn)化過程的空間關(guān)系,了解滲透系數(shù)非均質(zhì)性對(duì)潛流交換的影響、溶質(zhì)在空間上遷移轉(zhuǎn)化的差異和硝酸鹽去除率的變化。

1 非均質(zhì)潛流帶模型的建立

1.1 概念模型建立

使用COMSOL Multiphysics構(gòu)建二維潛流帶概念模型,用于模擬河流波動(dòng)引起的潛流帶內(nèi)氮素遷移轉(zhuǎn)化過程。模型幾何形狀取自Hornsby Bend站點(diǎn)上用全站儀測(cè)量的從河漫灘到河道中心的河岸高程輪廓,具體見圖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])

1.2 數(shù)學(xué)模型建立

本研究采用Richard方程和對(duì)流-彌散-反應(yīng)方程刻畫潛流帶內(nèi)水流和溶質(zhì)的遷移過程,其控制方程如下:

(1)

(2)

式中:p為壓力(Pa);g為重力加速度(m/s2);ρ為流體密度(kg/m3);μ為流體動(dòng)力黏度(Pa·s);Cm為比水容量(1/m);Se為有效飽和度(無量綱);S為儲(chǔ)水系數(shù)(1/Pa);ks為飽和水力滲透率(m2);kr為相對(duì)滲透率(無量綱),是有效飽和度的函數(shù);t為時(shí)間(s);z為位置高度(m);Qm為源匯項(xiàng)[kg/(m3·s)];θ為含水率(m3/m3);Ci為物質(zhì)i的濃度(mg/L);u為達(dá)西流速(m/s);D為水動(dòng)力彌散系數(shù)(m2/s);Ri為物質(zhì)i的反應(yīng)速率[mg/(L·d)]。

有氧呼吸作用 CH2O+O2→CO2+H2O

原位有機(jī)物或總有機(jī)碳的溶解可能導(dǎo)致DOC的輸入,該過程通過一級(jí)傳質(zhì)過程來模擬[17,37-40]:

(3)

式中:CDOC和CTOC分別為沉積物中DOC和總有機(jī)碳(TOC)的濃度(mg/L);α為一階傳質(zhì)系數(shù)(1/h);Kd為DOC分布系數(shù)(L/kg)。

有氧呼吸作用、硝化作用和反硝化作用由多重Monod動(dòng)力學(xué)描述[41]。存在DO時(shí),使用非競(jìng)爭(zhēng)性來抑制反硝化作用[42]。此外,系統(tǒng)中的氧通常用于能量?jī)?yōu)先級(jí)更高的反應(yīng),如有氧呼吸;而剩余的氧則作為具有較低能量?jī)?yōu)先級(jí)的二級(jí)溶質(zhì)被消耗,如反硝化反應(yīng)。而溶解氧需求的分布系數(shù)yO2是基于氧化和硝化反應(yīng)的吉布斯自由能計(jì)算的。綜上所述,各化學(xué)反應(yīng)項(xiàng)表示如下:

(4)

(5)

(6)

(7)

1.3 模型網(wǎng)格剖分與邊界設(shè)置

基于概念模型和數(shù)學(xué)模型,首先將潛流帶模型離散化。為了提高模型模擬精度和減少計(jì)算量,采用三角形網(wǎng)格進(jìn)行剖分(圖2),上邊界附近最小空間步長(zhǎng)為0.1 m,遠(yuǎn)離上邊界最大空間步長(zhǎng)為0.3 m,將模型離散為10 184個(gè)網(wǎng)格。模型模擬總時(shí)長(zhǎng)為24 h,時(shí)間步長(zhǎng)為1 h。

對(duì)于水流模型,河道中心剖面與底部邊界處理為隔水邊界;右側(cè)邊界和上邊界分別設(shè)置為定水頭邊界和透水層邊界,不考慮降雨或蒸發(fā)情況(具體設(shè)置參考Cardenas等[43])。河水位波動(dòng)近似表示成振幅為0.5 m、周期為24 h的正弦函數(shù)。

對(duì)于溶質(zhì)運(yùn)移模型,河道中心剖面與底部邊界處理為零通量邊界;右側(cè)邊界設(shè)置為流出邊界;上邊界設(shè)置為開放邊界,在這里允許發(fā)生溶質(zhì)的流入和流出,且外部濃度等于河流中溶質(zhì)濃度。假定在地下水水位上升時(shí),與大氣相連的非飽和帶中的DO處于飽和狀態(tài)。

表1 水流場(chǎng)和溶質(zhì)運(yùn)移場(chǎng)主要參數(shù)設(shè)置

2 結(jié)果與討論

2.1 潛流交換通量和流場(chǎng)分布情況

滲透系數(shù)的非均質(zhì)性影響河流-潛流帶的水動(dòng)力交換和流場(chǎng)分布,一個(gè)24 h周期內(nèi)河水位波動(dòng)和潛流交換通量,見圖3。

注:正值表示含水層向河流補(bǔ)給(出滲),負(fù)值表示河水向含水層補(bǔ)給(入滲);平均河水位和零通量則分別用紅色和藍(lán)色虛線表示。圖3 一個(gè)24 h周期內(nèi)河水位波動(dòng)(藍(lán)色實(shí)線)和潛流 交換通量(點(diǎn)線)Fig.3 Fluctuations in river level (solid blue line) and hyporheic exchange flux (dotted line) in a 24-hour period

由圖3可知:河水位上升時(shí),河水向含水層入滲(負(fù)值),河水位下降至平均河水位以下后,水由河水入滲含水層(負(fù)值)轉(zhuǎn)變?yōu)榱鞒龊畬?正值)向河水補(bǔ)給,并在18 h時(shí)河水位達(dá)到最低;18 h后,河水位雖然有所上升,但此時(shí)河水位仍低于平均河水位,含水層入滲補(bǔ)給河水;河流-含水層交界面上的入滲和出滲通量隨著河水位波動(dòng)而上升和下降,在河段河水位達(dá)到最大值前入滲流量達(dá)到峰值,而出滲通量在河段河水位達(dá)到最小值前達(dá)到峰值,這一趨勢(shì)與 Musial等[9]受潮汐影響的近溪流實(shí)地觀察結(jié)果和Shuai等[17]的模擬結(jié)果一致。

注:黑色箭頭為流線。圖4 不同方差下潛流帶滲透系數(shù)場(chǎng)和流場(chǎng)分布圖Fig.4 Diagram of permeability coefficient fields and flow fields of hyporheic zones under different variance values

由圖4可知,隨著潛流帶含水層非均質(zhì)性的增強(qiáng),流線越曲折,此外,含水層的非均質(zhì)性還影響潛流帶的交換速率,含水層的非均質(zhì)性越強(qiáng),潛流帶交界面處的等水位線向兩側(cè)延申越遠(yuǎn),潛流帶交換范圍越廣,反應(yīng)越劇烈。

2.2 溶質(zhì)濃度和生物地球化學(xué)反應(yīng)速率

圖5 不同方差下潛流帶中溶質(zhì)運(yùn)移過程圖Fig.5 Solute transport process in hyporheic zones with different variance values

圖6 不同方差下潛流帶中觀測(cè)點(diǎn)處溶質(zhì)濃度和化學(xué)反應(yīng)速率對(duì)比圖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ì)向潛流帶深處運(yùn)移,這在一定程度上了改變了硝化和反硝化帶的形態(tài),擴(kuò)大了反應(yīng)的作用區(qū)。與潛流帶內(nèi)溶質(zhì)濃度的空間分布形狀相類似,隨著含水層非均質(zhì)性的增強(qiáng),潛流帶內(nèi)高滲透層增加,硝化作用和反硝化作用區(qū)分布形狀越不規(guī)則,邊緣越“曲折”。

2.3 硝酸鹽去除率

為了定量分析潛流帶對(duì)硝酸鹽的去除情況,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ì)滲透深度均有所增加。總得來說,一個(gè)周期內(nèi)河流向潛流帶輸入的硝酸鹽總量和經(jīng)反硝化作用去除的硝酸鹽總量也隨之增加,潛流帶內(nèi)硝酸鹽去除率隨著滲透系數(shù)非均質(zhì)性程度的增加而增大。

3 結(jié) 論

本研究利用數(shù)值模擬手段研究了滲透系數(shù)非均質(zhì)性對(duì)潛流帶氮素遷移轉(zhuǎn)化和硝酸鹽去除的影響。主要得到結(jié)論如下:

1) 潛流內(nèi)滲透系數(shù)非均質(zhì)性程度的增加加快了潛流交換并提供了更高的反應(yīng)速率。隨著非均質(zhì)性程度的增加,潛流帶滲流面上入滲率和出滲率也越大,河水-含水層間的水力交換越劇烈。

2) 潛流內(nèi)滲透系數(shù)非均質(zhì)性影響了溶質(zhì)在含水層中的運(yùn)移情況,控制了潛流帶內(nèi)的反應(yīng)動(dòng)力學(xué),改變了溶質(zhì)分布形狀和硝化、反硝化作用的發(fā)生區(qū)域,隨著非均質(zhì)性的加強(qiáng),溶質(zhì)濃度和反應(yīng)區(qū)在空間的分布越不規(guī)則,分布區(qū)邊緣越“曲折”。

3) 滲透系數(shù)的非均質(zhì)性影響潛流帶去除硝酸鹽的能力,更高的滲透系數(shù)非均質(zhì)性表現(xiàn)出更好的硝酸鹽去除效果,硝酸鹽去除率也有所增加。