李曉宇,任仲宇,李芳春,井 琦*,魏明海
(1.北京工業(yè)大學建筑工程學院,北京 100124;2.水質科學與水環(huán)境恢復工程北京市重點實驗室,北京 100124;3.同方知網(北京)技術有限公司,北京 100192;4.生態(tài)環(huán)境部環(huán)境規(guī)劃院,北京 100012)
阿特拉津(Atrazine)是目前世界上廣泛使用的除草劑之一,也是土壤和地下水中最常檢測到的農藥之一[1],已在中國的松花江流域、遼河流域和長江流域中檢測到[2-4]。阿特拉津被認定為有潛在的致癌作用并可影響人體的內分泌系統(tǒng),美國、歐共體和日本等均把它列入內分泌干擾物名單[5-6]。阿特拉津及其代謝產物可在水和土壤中存在數(shù)十年[7],其殘留會對飲用水源和生態(tài)環(huán)境構成威脅,或通過食物鏈傳遞和生物富集作用對人體健康造成危害。因此,研究阿特拉津在多孔介質中的遷移與變化規(guī)律具有十分重要的現(xiàn)實意義。
國內外學者對影響阿特拉津遷移的主要過程(吸附與解吸)及作用機理的相關研究已有大量報道[8-9]。Kovaios等[10]研究了土壤對阿特拉津的吸附作用,發(fā)現(xiàn)阿特拉津的吸附是可逆的,并觀察了吸附和解吸過程。針對土壤中阿特拉津的易遷移性,許多學者進行了室內土柱試驗和田間試驗,并在此基礎上通過模型模擬阿特拉津的遷移,從而評價其對地下水的污染風險[4]。Mao等[11]、Zheng等[12]利用兩點(兩區(qū))非平衡吸附模型模擬阿特拉津在土壤中的遷移,并取得了較好的模擬效果。Prata等[13]采用三點非平衡吸附模型成功地擬合了淋溶條件下阿特拉津在巴西一種氧化土中的穿透曲線。毛萌等[14]建立了包含吸附項和降解項的對流彌散方程,模擬室內滴灌施藥條件下阿特拉津在土壤中的遷移。Van Genuchten等[15]發(fā)展了兩點模型的控制方程,在原方程基礎上增加了降解項。Gamerdinger等[16]應用帶有降解項的兩點模型模擬了室內定流速條件下阿特拉津在土柱中的遷移。
綜上所述,已有研究對模型的改進主要表現(xiàn)在改進模型源匯項(例如考慮吸附、降解項等)或將兩點非平衡吸附模型改進為三點非平衡吸附模型。單點非平衡吸附模型是兩點非平衡吸附模型的特例[15],Selim等[17]曾利用單點非平衡吸附模型模擬分析2,4-D除草劑在飽和及非飽和分層土壤中的遷移,但未討論該模型的模擬效果。吸附-解吸模型可用來模擬病毒、細菌和膠體等在多孔介質中的遷移[18],Yu等[19]還應用吸附-解吸模型模擬了納米零價鐵在飽和砂柱中的遷移,杜青青等[20]應用吸附-解吸模型模擬了氨氮在污染場地的遷移,但有關吸附-解吸模型模擬分析阿特拉津等農藥在多孔介質中的遷移報道甚少。筆者在利用單點非平衡吸附模型模擬阿特拉津遷移過程中,發(fā)現(xiàn)在淋洗階段誤差較大,進而考慮采用吸附-解吸模型對阿特拉津遷移進行模擬,通過改進吸附、解吸的源匯項表達,以期消減此誤差。
本文通過構建室內土柱阿特拉津動態(tài)遷移試驗,獲取模擬所需參數(shù),利用HYDRUS-1D軟件,分別建立了基于單點非平衡吸附模型和吸附-解吸模型的運移模型,采用整體調參和分段調參兩種模擬方式,模擬阿特拉津在飽和壤質砂土中的遷移與分布規(guī)律,深入分析了兩種模型模擬結果差異背后的物理機理,并對源匯項的分項逐項進行計算,定量分析了兩種吸附模型對模擬結果的影響過程和影響程度,可為深刻認識阿特拉津的遷移特征及其環(huán)境風險評價提供科學依據(jù)。
1.1.1 材料與儀器
試驗材料:阿特拉津(純度98%,國家標準物質中心)、氯化鉀(分析純,北京化工廠)、甲醇(色譜純、國藥集團化學試劑有限公司)。供試土壤采自北京市永定河河灘,采樣深度40 cm,土壤取回后自然風干,過2 mm篩備用。用簡易比重計法進行了該土壤的粒徑分析,并按國際制將土壤質地定名為壤質砂土。分析土壤的砂粒含量為86.653%,粉粒含量為11.795%,黏粒含量為1.553%,有機質含量為4.230%,干容重為1.656 g·cm-3,有效孔隙度為0.130,pH值為6.300。
試驗儀器:電子天平(JA2003,上海恒平科學儀器有限公司)、電導率儀(sensION5,哈希水質分析儀器有限公司)、高效液相色譜儀(Agilent1100,安捷倫科技有限公司)等。
模擬裝置如圖1所示,由進水裝置、土柱主體、出水裝置三部分組成。進水裝置由馬氏瓶、進水管和升降支架組成,馬氏瓶用于維持土柱定水頭淹水狀態(tài),升降支架用于承載馬氏瓶,可調節(jié)高度。土柱主體由有機玻璃柱和充填介質組成,有機玻璃柱工作長度36 cm,內徑6 cm,柱頂設置溢流孔,柱中及柱底分別設置取樣孔。出水裝置由出水管和集液箱組成。
圖1 柱試驗裝置示意圖Figure 1 Schematic diagramof column experimental device
1.1.2 土柱遷移試驗
在土柱底部鋪設5 cm厚粒徑為0.5~2 cm的卵礫石,然后將風干過篩后的土壤按10 cm層厚分層填入土柱,裝填過程中對土壤層間打毛并在表層土壤上部鋪設3 cm厚粒徑為0.5~2 cm的卵礫石。土柱填裝完成后對土柱進行排氣飽和并調整進水出水的水位差,得到穩(wěn)定滲流速度,在維持水位差恒定條件下將濃度為22μg·mL-1的阿特拉津滲濾液從頂端注入柱內,注入時長200 min,總注入水量1650 cm3,約12.5個土壤有效孔隙體積(12.5 PV),注入結束后瞬時用去離子水淋洗土柱,淋洗時間150 min,總淋洗水量1238 cm3(9.4 PV)。定時從取樣孔取樣,并測定滲流液中阿特拉津濃度。每次每個觀測孔取樣2個,測定結果取樣品的平均值。
1.1.3 測定方法
應用安捷倫高效液相色譜儀測定阿特拉津濃度。試驗采用C18色譜柱,設置流動相(甲醇∶水=70∶30)的流速為1 mL·min-1,在225 nm波長條件下用UVVIS檢測器測定。通過比對峰面積/濃度標線,計算得到不同時刻阿特拉津濃度。
HYDRUS軟件可用于模擬分析非飽和、飽和多孔介質中的水和溶質運移[18]。本文利用HYDRUS-1D軟件建立土壤水流和溶質運移數(shù)值模型,模擬阿特拉津溶液入滲進入飽和土柱后的遷移過程以及隨后淋洗條件下阿特拉津在飽和土柱中的遷移過程,設定時間單位為min,質量單位為μg,長度單位為cm。根據(jù)土柱試驗的運行情況,將模擬時長定為350 min。
1.2.1 水流模型的基本方程
一維平衡水流運動采用Richards方程來描述[18,21]:
式中:h為土壤壓力水頭,cm;θ為土壤體積含水率;t為時間,min;S為源匯項,min-1;α為水流方向與垂直方向夾角(即:α=0°代表垂向流動,90°代表水平流動,0°<α<90°代表傾斜流);x為空間位置,cm;K為非飽和水力傳導度函數(shù),cm·min-1,如式(2):
式中:Kr為相對水力傳導度;Ks為飽和水力傳導度,cm·min-1。
本研究水流模型概化為均質各項同性飽和一維垂向穩(wěn)定流,因此方程(1)可簡化為下式:
式中各參數(shù)意義同上。
1.2.2 水流模型的初始條件及邊界條件
模型初始條件為整體飽和狀態(tài),即初始含水率為飽和含水率。
模型上邊界為供水邊界,下邊界為排水邊界,且均為定水頭邊界,水頭差為5.17 cm,通過插值得到土柱內的初始流場。
1.2.3 溶質遷移模型的基本方程
通過改進的對流彌散方程來模擬阿特拉津在飽和多孔介質中的遷移,該方程[22]用一個源匯項表示阿特拉津與固體基質的反應(吸附和解吸過程)。
式中:n為孔隙度;C為溶質濃度,μg·mL-1;S為多孔介質中溶質固相濃度,μg·g-1;D為彌散系數(shù),cm2·min-1;v為達西速度,cm·min-1;x為空間位置,cm;ρ為多孔介質干容重,g·cm-3;t為時間,min。
根據(jù)不同的吸附機理,源匯項方程式有不同的表達形式。本文分別采用單點非平衡吸附模型和吸附-解吸模型對比模擬分析阿特拉津在土柱中的遷移過程。
式中:α是一級反應物質交換率系數(shù)(表示單位時間液相與固相間溶質質量轉化大小),min-1;Kd是Henry線性吸附系數(shù),mL·g-1。
式中:ka為吸附系數(shù),min-1;kb為解吸系數(shù),min-1;ψ 為溶質吸附函數(shù);Smax是固相最大吸附量,μg·g-1;θ為土壤體積含水率,在飽和狀態(tài)下,砂土中θ可近似認為等于孔隙度。
1.2.4 溶質遷移模型的初始條件及邊界條件
模型初始條件為初始時刻土柱各位置溶質濃度為0。
模型上邊界為濃度通量邊界,下邊界設定0梯度邊界條件。佩克萊數(shù)Pe表示對流與擴散在溶質遷移中所占比例[24]。Parlange等[25]認為,當Pe>4時,近似認為多孔介質流體中溶質濃度與某時段內按出流量計算的平均溶質濃度相等,此時邊界條件為0濃度梯度邊界條件。Pe計算公式如下:
式中:L為柱長,cm;u為實際平均速度,u=v/n,cm·min-1;D為彌散系數(shù),cm2·min-1。通過計算本試驗佩克萊數(shù)Pe=49.860>4,符合Parlange等[25]研究的判定條件,即擴散作用相對于對流作用在溶質遷移中所占的比例較小,此時出流溶質的濃度梯度近似為0,故下邊界設定為零梯度邊界。
本模擬中需要的參數(shù)有土壤水力特征參數(shù)和溶質特征參數(shù),土壤特征參數(shù)主要指孔隙度、飽和滲透系數(shù)等;溶質特征參數(shù)主要指彌散系數(shù)、吸附系數(shù)等。
1.3.1 土壤水力特征參數(shù)
土壤水力特性參數(shù)的準確性對于反映實際土壤水分運動過程具有重要意義[26]。本研究采用常水頭法測定土壤介質飽和滲透系數(shù)Ks,容積排水法測定土壤介質孔隙度。測量結果土壤飽和滲透系數(shù)為0.340 cm·min-1,孔隙度為0.349。
1.3.2 彌散系數(shù)
彌散系數(shù)是表征可溶性物質通過滲透介質時彌散現(xiàn)象強弱的指標,可通過實驗利用投加示蹤劑測定其穿透曲線的方式獲取。本研究在室內土柱模擬裝置中,以Cl-作為示蹤劑,在飽和一維穩(wěn)定流條件下淋濾KCl溶液,并在土柱出口端定時測定Cl-濃度,利用公式(8)[27]計算彌散系數(shù)。計算得到彌散系數(shù)為0.600 cm2·min-1。
彌散系數(shù)計算表達式:
式中:t0.16、t0.5、t0.84為柱中一個固定點x處相對濃度,分別為0.16、0.5、0.84時對應的時間,min。
Cl-濃度測試采用電導率法。配制不同濃度的KCl溶液,并測定其對應的電導率值,繪制電導率與KCl濃度的關系曲線,試驗時用電導率儀測定樣品電導率,利用電導率與KCl濃度關系將電導率換算為濃度。
1.3.3 等溫吸附系數(shù)
吸附系數(shù)的確定由室內吸附試驗獲得,本研究中阿特拉津在土壤中的吸附行為符合Henry等溫線性吸附,線性吸附系數(shù)Kd為0.372 mL·g-1;將土壤在濃度為22μg·mL-1的阿特拉津溶液中恒溫振蕩200 min,得到吸附量為6.450μg·g-1,以此作為本試驗條件下的固相最大吸附量Smax。
在水力試驗基礎上確定水流模型參數(shù),建立土柱一維水流穩(wěn)定流模型,調整邊界水頭以改變土柱內的流場,將觀測孔的實測測壓管水頭值和計算水頭值進行對比,反復進行三次,取誤差較小的情形所對應的參數(shù)為模型水流參數(shù)。模型中涉及的土壤水力特征參數(shù)如飽和滲透系數(shù)為0.340 cm·min-1,孔隙度為0.349。
水質模型中源匯項涉及的參數(shù)有彌散系數(shù)D、一級反應物質交換率系數(shù)α、線性吸附常數(shù)kd、吸附系數(shù)ka、解吸系數(shù)kb、固相最大吸附量Smax。彌散系數(shù)D、線性吸附常數(shù)kd、固相最大吸附量Smax通過實驗得到,并在模型中適當調整,最終得到這些參數(shù)值:D為0.600 cm2·min-1,kd為0.372 mL·g-1,Smax為6.450 μg·g-1;針對以上兩種吸附模型中不易通過試驗直接測得的參數(shù)如一級反應物質交換率系數(shù)α、吸附系數(shù)ka、解吸系數(shù)kb進行調參,得到綜合反映整個試驗過程的“平均的”α、ka和kb,分別為 0.200、0.179 min-1和 0.126×10-5min-1。
實際情況下,隨注入、淋洗階段時間的增長,模型中的一級反應物質交換率系數(shù)、吸附和解吸系數(shù)會發(fā)生變化。為了探究實際情況下各反應系數(shù)的變化過程及對模擬效果的影響,本研究以觀測孔1和觀測孔2為調參檢驗點,結合阿特拉津穿透曲線形態(tài),將模擬按照穿透曲線的“上升段”“峰值段”“下降段”分段進行調參模擬,得到各階段參數(shù)值如表1所示。
圖2、圖3分別為取樣孔1和取樣孔2阿特拉津濃度實測值與整體調參(未分段調參)擬合的單點非平衡吸附模型、吸附-解吸模型模擬計算結果對比圖。圖4、圖5分別為取樣孔1和取樣孔2阿特拉津濃度實測值與分段調參擬合的單點非平衡吸附模型、吸附-解吸模型模擬計算結果對比圖。
單點非平衡吸附模型在阿特拉津注入過程中擬合情況較好,各觀測點濃度模擬值和實測值較為一致,但在淋洗過程中,模擬結果相比于試驗結果有較為明顯的誤差,不能較好地重現(xiàn)整個過程(見圖2和圖3);吸附-解吸模型在整個過程中模擬值與實測值趨勢更為一致,且誤差較小,可以較好地描述阿特拉津在土壤中的遷移過程。通過觀察圖4、圖5并與圖2、圖3對比可知,分段調參擬合的單點非平衡吸附模型的擬合度在注入階段提升不明顯,在淋洗階段有顯著提升;分段調參擬合的吸附-解吸模型的擬合度在注入階段和淋洗階段均提升很小。
選用Nash-Suttcliffe模擬效果系數(shù)(NSC)和均方根誤差(RMSE)評價模型計算結果,計算公式分別為:
表1 各階段一級反應物質交換率系數(shù)、吸附和解吸系數(shù)擬合值Table 1 Fitting values of exchange rate coefficient of first order reactive substance,adsorption and desorption coefficients at each stage
圖2 觀測孔1整體調參擬合:單點非平衡吸附模型與吸附-解吸模型模擬結果對比圖Figure 2 Integral parameter fitting of observation hole 1:The simulation result compartion chart of one-site non-equilibriumadsorption model and adsorption-desorption model
圖3 觀測孔2整體調參擬合:單點非平衡吸附模型與吸附-解吸模型模擬結果對比圖Figure 3 Integral parameter fitting of observation hole 2:The simulation result compartion chart of one-site non-equilibriumadsorption model and adsorption-desorption model
圖4 觀測孔1分段調參擬合:單點非平衡吸附模型與吸附-解吸模型模擬結果對比圖Figure 4 Segmental parameter fitting of observation hole 1:The simulation result compartion chart of one-site non-equilibriumadsorption model and adsorption-desorption model
式中:Xobs為實測值,Xmodel為模型計算值,Xˉobs為實測值的算術平均值。當模型計算值和實際監(jiān)測值相等時,NSC=1,模擬效果最好;通常NSC在0~1之間,NSC越大,說明計算值與觀測值匹配程度越好,一般當NSC>0.75時,表示模擬結果可以接受。
式中:N表示實測值個數(shù),Xobs,i和Xmodel,i為試驗觀測值與模型計算值。RMSE越小,說明計算值與觀測值匹配程度越好。
分整體調參擬合和分段調參擬合兩種情形,按照阿特拉津注入階段、淋洗階段和注入-淋洗全階段,計算單點非平衡吸附模型和吸附-解吸模型的NSC和RMSE值,如表2所示。
由表2可知,整體調參擬合情形下,單點非平衡吸附模型模擬注入階段(吸附階段)NSCs與RMSEs分別為0.909~0.922和2.752~3.167,吻合較好;而對于淋洗階段,NSCs為負,說明擬合度很低;單點非平衡吸附模型模擬阿特拉津注入-淋洗的全過程的NSCs與RMSEs分別為0.138~0.593和6.180~8.395。分段調參擬合后,注入階段的NSCs提高到0.919~0.941,RMSEs減小到2.405~2.986,淋洗階段的NSCs提高到0.927~0.940,RMSEs減小到2.309~2.036,注入-淋洗全階段的NSCs提高到 0.921~0.941,RMSEs減小到 2.204~2.717。說明用分段擬合的單相非平衡吸附模型模擬阿特拉津注入、淋洗和注入-淋洗的全過程均是適宜的。
同樣通過觀察吸附-解吸模型的孔1和孔2的NSC和RMSE值可知,無論整體調參模擬和分段調參模擬,吸附-解吸模型模擬阿特拉津注入、淋洗和注入-淋洗全過程均是較為適宜的。分段擬合的吸附-解吸模型模擬阿特拉津注入、注入-淋洗全過程的模擬效果系數(shù)有微弱提高。
圖5 觀測孔2分段調參擬合:單點非平衡吸附模型與吸附-解吸模型模擬結果對比圖Figure 5 Segmental parameter fitting of observation hole 2:The simulation result compartion chart of one-site non-equilibrium adsorption model and adsorption-desorption model
單點非平衡吸附模型認為固相中的吸附量S與一級反應物質交換率系數(shù)α和線性吸附常數(shù)與溶質液相濃度的乘積KdC有關。S和KdC隨著吸附過程的發(fā)生不斷發(fā)生變化。根據(jù)單點非平衡吸附模型,可將阿特拉津在土柱中的遷移、吸附過程描述如下:初始時刻柱中各位置阿特拉津固相及液相濃度均為0,阿特拉津隨水流遷移進入土柱,土柱中液相濃度升高,同時發(fā)生吸附作用,固相濃度隨之增加,達到吸附平衡后停止吸附。隨后,注入去離子水,吸附在固相中的阿特拉津逐漸解吸到水中隨水流出。由圖2、圖3可知,阿特拉津連續(xù)注入階段(吸附階段)模擬值與實測值偏差較小,但在去離子水注入階段(淋洗階段)則偏差很大。這是因為整個注入-淋洗階段的一級反應物質交換率系數(shù)α為同一常數(shù),淋洗階段由于采用了與注入階段相同的α,導致淋洗階段實測值與模擬值差異較大。分段擬合后,單點非平衡吸附模型由于在注入階段和淋洗階段采用了與對應階段相匹配的α,因此注入階段和淋洗階段的擬合度均較高。
表2 各階段模型模擬效果系數(shù)(NSC)和均方根誤差(RMSE)值Table 2 Simulation effect coefficient(NSC)and root-mean-squared error(RMSE)values of models at each stage
為了進一步說明兩種吸附模型的結果差異,本文對兩種吸附模型在不同階段的源匯項值進行了準確計算,分別列示于表3與表4。
表3列出了觀測孔1單點非平衡吸附模型源匯項計算值。由表3可知:隨著阿特拉津溶液的持續(xù)注入,模型計算得到的液相濃度和固相濃度逐漸增大,均為正值且逐漸減小,到阿特拉津注為0.003 μg·cm-3·min-1,此時固相濃度仍有較小的正向增長;在淋洗階段,隨著去離子水的持續(xù)注入,模型計算得到的液相濃度和固相濃變?yōu)樨撝担ㄘ撝当硎窘馕^程)且絕對值與注入阿特拉津時的源匯項值相近,即吸附到固相中的阿特拉津逐漸解吸到溶液中,解吸速率與注入阿特拉津溶液時的吸附速率大致相當,因此采用單點非平衡吸附模型在淋洗階段模擬值遠超出實測值。
吸附-解吸模型認為吸附和解吸速率不同,將模型分為吸附項和解吸項分別描述,用ka和kb代表吸附和解吸的速率。在吸附項中,吸附量與溶質液相濃度C和最大固相濃度Smax有關,當吸附量達到Smax時停止吸附。根據(jù)吸附-解吸模型,可將阿特拉津在土柱中遷移、吸附過程描述如下:初始時刻柱中各位置阿特拉津液相和固相濃度均為0,隨著阿特拉津注入土柱,土柱中液相濃度升高,同時發(fā)生吸附作用,固相濃度隨之增加,隨著土柱中各位置吸附量達到最大值后,吸附不再發(fā)生,整個過程中阿特拉津也會發(fā)生少量的解吸;對于淋洗階段,同樣存在著吸附和解吸作用,解吸系數(shù)kb在數(shù)值上很小,說明吸附到壤質砂土中的阿特拉津解吸速率很慢,這也是吸附-解吸模型計算值遠低于單點非平衡模型計算值的原因。
表4列出了觀測孔1吸附-解吸模型源匯項計算值。由表4可知:隨著阿特拉津溶液的持續(xù)注入,模型計算的液相濃度和固相濃度逐漸增大,分別達到注入濃度和固相最大吸附量后維持不變,源匯項隨著固相濃度的持續(xù)增大逐漸減小,達到最大固相吸附量后變成負值且維持恒定,此時的源匯項為-1.388×10-5μg·cm-3·min-1,可知固相濃度雖然在逐漸降低,但下降速率近似為0;在淋洗階段,隨著去離子水的持續(xù)注入,模型計算的液相濃度逐漸降低,源匯項為-1.388×10-5μg·cm-3·min-1,表明吸附在固相上的阿特拉津解吸速率極慢。通過觀察阿特拉津液相濃度實測穿透曲線發(fā)現(xiàn)沒有拖尾現(xiàn)象,也說明了其實際解吸速率很小。對比表3、表4的源匯項值計算結果,在解吸階段單點非平衡吸附模型的源匯項計算值比吸附-解吸模型的源匯項計算值高出4個數(shù)量級,這是引起兩種模型計算結果差異較大的原因。
表3 觀測孔1單點非平衡吸附模型源匯項計算值Table 3 Calculation values of source and sink terms of the onesite non-equilibriumadsorption model of the observation hole 1
表4 觀測孔1吸附-解吸模型源匯項計算值Table 4 Calculation values of source and sink terms of the adsorption-desorption model of the observation hole 1
(1)整體調參下,單點非平衡吸附模型在連續(xù)注入阿特拉津的條件下,飽和壤質砂土中阿特拉津遷移模擬的NSCs系數(shù)為 0.909~0.922,RMSEs為 2.752~3.167,采用單點非平衡吸附模型是適合的;淋洗條件下阿特拉津遷移模擬的NSCs為負值,說明單點非平衡吸附模型不適用于淋洗條件下阿特拉津遷移模擬。
(2)整體調參下,吸附-解吸模型模擬連續(xù)注入、淋洗過程中阿特拉津在飽和壤質砂土中的遷移模擬的NSCs為 0.901~0.954,RMSEs為 2.037~3.289,表明吸附-解吸模型模擬連續(xù)注入、淋洗過程中阿特拉津遷移均是適合的。
(3)分段調參可提高單點非平衡吸附模型模擬阿特拉津連續(xù)注入、淋洗過程中在飽和壤質砂土中遷移的擬合度。