HYDRUS-1D軟件在土壤環(huán)境影響評(píng)價(jià)中的應(yīng)用

丁素玲

(中材地質(zhì)工程勘查研究院有限公司，北京 100102)

為加強(qiáng)土壤污染防治、改善土壤環(huán)境質(zhì)量、保障農(nóng)產(chǎn)品質(zhì)量和人居環(huán)境安全，國(guó)家和地方陸續(xù)出臺(tái)了各種與土壤環(huán)境保護(hù)相關(guān)的法律法規(guī)，并于2018年9月由生態(tài)環(huán)境部發(fā)布了《環(huán)境影響評(píng)價(jià)技術(shù)導(dǎo)則 土壤環(huán)境(試行)》(HJ 964—2018)(以下簡(jiǎn)稱《導(dǎo)則》)，該《導(dǎo)則》規(guī)定了土壤環(huán)境影響評(píng)價(jià)的內(nèi)容和方法，同時(shí)要求選擇適宜的預(yù)測(cè)方法，預(yù)測(cè)評(píng)價(jià)建設(shè)項(xiàng)目各實(shí)施階段不同環(huán)節(jié)與不同環(huán)境影響防控措施下的土壤環(huán)境影響，但《導(dǎo)則》中未給出具體的土壤環(huán)境影響預(yù)測(cè)方法。本文以某飛灰固化填埋項(xiàng)目為例，在土壤理化性質(zhì)分析[1]和工程試驗(yàn)的基礎(chǔ)上，借助HYDRUS-1D軟件建立數(shù)值模型，模擬污染物在土壤中的運(yùn)移過程，為土壤環(huán)境影響預(yù)測(cè)提供一種簡(jiǎn)便可行的方法。

1 研究區(qū)土壤及污染源概況

研究區(qū)土壤主要為磚紅壤、黏土和強(qiáng)風(fēng)化玄武巖，團(tuán)粒結(jié)構(gòu)，含砂量約為20%，總厚度為9.0 m。其中，磚紅壤厚度為1.5 m，埋深0 m～1.5 m；黏土層厚度為0.5 m，埋深1.5 m～2.0 m；強(qiáng)風(fēng)化玄武巖厚度為7.0 m，埋深2.0 m～9.0 m。

研究對(duì)象的污染源為飛灰固化物填埋場(chǎng)滲濾液收集池。飛灰是指生活垃圾焚燒煙氣凈化系統(tǒng)捕集物和煙道及煙囪底部沉積的底灰，其成分比較復(fù)雜，變化范圍也較大。其中大部分是硅酸鹽和含鈣、鋁、鐵、鉀等金屬的化學(xué)物質(zhì)，另外還含有多種重金屬如Hg、Pb、Cd、As、Zn、Sb、 Se等，以及其他有毒有機(jī)物如二英等。飛灰填埋場(chǎng)滲濾液主要來源為降雨，其主要成分為pH、CODcr、BOD5、氨氮、總磷、砷、總鉻、六價(jià)鉻、鉛、汞、鎳、銅、鎘等重金屬。經(jīng)對(duì)滲濾液中各污染物濃度的分析和比較可知：本次研究對(duì)象飛灰固化物填埋場(chǎng)的滲濾液中污染物砷的指標(biāo)含量最高，因此，將砷作為預(yù)測(cè)因子，其濃度為0.224 mg/L。污染物對(duì)土壤造成污染的途徑主要為經(jīng)地表徑流或下滲進(jìn)入土壤[2-3]。

2 土壤中水分和溶質(zhì)運(yùn)移原理

水是污染物運(yùn)移的載體，若進(jìn)行土壤中污染物運(yùn)移的模擬需要先對(duì)土壤中水分運(yùn)移過程進(jìn)行刻畫。

2.1 水流運(yùn)動(dòng)方程

一般認(rèn)為，水在土壤中運(yùn)移符合活塞流模式。土壤中水流運(yùn)動(dòng)的控制方程為一維垂向飽和——非飽和水分運(yùn)動(dòng)方程(Richards方程)：

(1)

式中，θ(h)為土壤體積含水量；h為壓力水頭[L]，飽和帶大于零，非飽和帶小于零；z、t分別為垂直方向坐標(biāo)變量[L]、時(shí)間變量[T]；k(h)為垂直方向的水力傳導(dǎo)度[LT-1]；s為作物根系吸水率[T-1]，不考慮根系吸水時(shí)，取s=0。

采用無滯后效應(yīng)的Van Genuchten-Mualem模型，它嵌入了Scott(1983)、Kool和Parker(1987)經(jīng)驗(yàn)?zāi)Ｐ椭械募俣ǎ何鼭?脫濕)掃描線與主吸濕(脫濕)曲線成比例變化，并運(yùn)用一個(gè)比例程序，將用戶定義的水力傳導(dǎo)曲線與參考土壤相比較，通過線性比例變換，獲得給定土壤剖面上的近似水力傳導(dǎo)變量。

方程(1)中相關(guān)參數(shù)可用公式(2)和公式(3)進(jìn)行求解：

(2)

(3)

式中，θr為土壤殘余含水量；θs為土壤飽和含水量；α為土壤持水參數(shù)[L-1]；m、n均為土壤持水指數(shù)；ks為土壤飽水滲透系數(shù)[LT-1]；l為有效孔隙度；se為貯水率[L-1]。

在解決實(shí)際問題時(shí)，必須加上描述具體問題的初始和邊界條件，組成某種狀態(tài)下表述土壤水分運(yùn)動(dòng)和溶質(zhì)運(yùn)移的數(shù)學(xué)模型。各條件設(shè)置如下:

初始條件：h(z,0)=h0

式中，z為地表至下邊界距離[L]；qs為污水下滲通量[LT-1]；h(z，t)為土壤壓力水頭[L]。

2.2 溶質(zhì)運(yùn)移方程

溶質(zhì)運(yùn)移方程建立在水流模型的基礎(chǔ)上，不考慮溶液密度的變化，且本著風(fēng)險(xiǎn)最大的原則，忽略污染物吸附、解析和自然衰減等物理、化學(xué)、生物反應(yīng)，只考慮對(duì)流、彌散作用，建立的土壤中一維垂向溶質(zhì)運(yùn)移方程如下：

(4)

式中，θ為土壤體積含水率；c為土壤液相中污染物的濃度；D為綜合彌散系數(shù)；q為體積流動(dòng)通量密度。

初始條件：假定土壤中污染物的初始濃度為零。

c(z,t)=0t=0,L≤z<0，L為土柱深度。

上邊界條件：連續(xù)點(diǎn)源污染(污染物以固定濃度c0連續(xù)注入)的情景下，上邊界為給定濃度的第一類Dirichlet邊界。

c(z,t)=c0t>0，z=0

下邊界條件：由于模擬選擇的下邊界為潛水面，污染物呈自由入滲狀態(tài)，邊界內(nèi)外的濃度相等，故而將其設(shè)定為不存在彌散通量的第二類Neumann零梯度邊界。

通過求解上述數(shù)學(xué)模型便可定量刻畫污染物在土壤中的不同時(shí)刻、不同位置處的污染物濃度。

3 土壤中污染物運(yùn)移模擬過程

HYDRUS是由美國(guó)國(guó)家鹽土改良中心(US Salinity Laboratory)于1991年研制的用于模擬變飽和多孔介質(zhì)中水分、能量、溶質(zhì)運(yùn)移的數(shù)值模型，已得到廣泛的認(rèn)可與應(yīng)用[4-5]。該軟件能夠較好地模擬水分與污染物在土壤中的分布和運(yùn)移規(guī)律，也可以與其他地下水、地表水模型相結(jié)合，從宏觀上分析水資源問題[6]。HYDRUS-1D是HYDRUS模型的一維版本，是一種常用的水和溶質(zhì)模擬軟件。近年來該軟件得到廣泛的應(yīng)用，如葉永紅等[7]使用該軟件預(yù)測(cè)了蘭州西固商業(yè)石油儲(chǔ)備庫泄漏對(duì)包氣帶的影響，郝芳華等[8]將其應(yīng)用于河套灌區(qū)典型區(qū)域的土壤氮元素遷移，楊洋等[9]利用HYDRUS-1D模擬了填埋場(chǎng)滲濾液中氨氮在包氣帶中的遷移轉(zhuǎn)化規(guī)律等。

用HYDRUS-1D模擬土壤中污染物運(yùn)移過程的方法如下：

(1)目標(biāo)土層剖分及觀測(cè)點(diǎn)布置

在HYDRUS-1D(Soil Profile-Graphical Editor)模塊中對(duì)包氣帶土層進(jìn)行設(shè)定。將整個(gè)包氣帶剖面劃分為900層，每層1 cm，總厚度為9.0 m。其中，0 m～1.5 m土層為壤土層，1.5 m～2.0 m為黏土層，2.0 m～9.0 m為強(qiáng)風(fēng)化玄武巖層。在潛水面處布置1個(gè)觀測(cè)點(diǎn)N1。

(2)參數(shù)設(shè)置

HYDRUS-1D軟件所需參數(shù)分為兩類：一類是土壤特征參數(shù)，如土壤的飽和含水量、殘余含水量和土壤水力曲線參數(shù)等水分參數(shù)；另一類是運(yùn)移參數(shù)，如滲透系數(shù)、彌散系數(shù)等。

本模型不考慮土壤對(duì)污染物的吸附和污染物的轉(zhuǎn)化。各參數(shù)除土層飽和滲透系數(shù)使用注水試驗(yàn)的實(shí)測(cè)值外，其余各參數(shù)均采用HYDRUS-1D軟件自帶的相應(yīng)土層的經(jīng)驗(yàn)參數(shù)值。模型中主要參數(shù)見表1。

表1 模型主要參數(shù)表

(3)邊界條件和初始條件

上邊界：上邊界為滲濾液收集池的底部，假定池底破損，污染物持續(xù)入滲，則上邊界為固定濃度邊界，濃度C0=0.224 mg/L。

下邊界：假設(shè)潛水面隨時(shí)間的變化幅度較小，可忽略不計(jì)，則下邊界可以設(shè)為壓力水頭為零的第一類Dirichlet邊界。

初始條件：本次模擬初始水頭設(shè)定h0=-100 cm(土壤處于干燥狀態(tài))。土壤中砷的初始濃度C0=0 mg/L。

4 模擬實(shí)驗(yàn)結(jié)果分析

本次模擬對(duì)象為滲濾液收集池，假設(shè)池底防滲層破損，滲濾液通過破損區(qū)域泄漏到土壤中。將污染源概化為穩(wěn)定的點(diǎn)源污染，污染物砷以固定濃度0.224 mg/L持續(xù)入滲。泄漏的時(shí)長(zhǎng)是25年(9125天，為本項(xiàng)目的服務(wù)年限)。預(yù)測(cè)結(jié)果見圖1、圖2。

由以上模擬結(jié)果可知：

在泄漏發(fā)生后，飛灰固化物滲濾液向下遷移形成垂向污染暈，在對(duì)流和彌散的作用下，污染暈鋒面處的濃度隨深度加深而逐漸減小。隨著滲濾液不斷下滲，污染物的累積速度超過土壤的自凈速度和容納能力，土壤中污染深度和污染物濃度均不斷增加。

滲濾液泄漏后第5天，污染暈最大下滲深度為5.7 m，土壤水中砷的最大濃度為0.00001 mg/L；泄漏后第50天，污染暈最大下滲深度為9.0 m，此時(shí)砷已經(jīng)穿透包氣帶土層到達(dá)潛水面處，土壤水中砷濃度最大值為0.000233 mg/L；隨著污染物的持續(xù)入滲，泄漏后第365天，土壤水中砷濃度不斷增加，最大值為0.0818 mg/L；泄漏后第9125天，土壤水中砷濃度達(dá)到峰值，為0.224 mg/L；此后，滲濾液停止泄漏，土壤水中砷濃度逐漸降低，至第10 220天，土壤水中砷濃度最大值為0.0362 mg/L。

(注：N1控制點(diǎn)位置：潛水面)圖1 控制點(diǎn)N1處砷濃度隨時(shí)間變化圖Fig.1 As concentration change with time at control point N1

(注：T1～T6分別代表泄漏開始后第1天、5天、50天、365天、9125天、10 220天)圖2 不同時(shí)間節(jié)點(diǎn)下砷濃度隨深度變化圖Fig.2 As concentration change with depth at different time nodes

從潛水面處N1控制點(diǎn)的污染物濃度變化預(yù)測(cè)結(jié)果可知：污染物泄漏前，N1處土壤水中污染物砷濃度為0 mg/L；在飛灰固化物滲濾液收集池泄漏后的第12.4天，N1處土壤水中砷濃度為0.00001 mg/L，這說明砷的污染鋒面已經(jīng)到達(dá)此處，污染物已穿過整個(gè)包氣帶土壤層到達(dá)地下水中；泄漏后的第125天，N1處土壤水中砷的濃度為0.01118 mg/L；泄漏后第3189天，N1處土壤水的砷濃度達(dá)到峰值且持續(xù)5970天，直至第9159天(污染源消失后的第34天)砷濃度逐漸下降。

5 結(jié)論與建議

本文以某飛灰固化物填埋項(xiàng)目為例，假設(shè)飛灰固化物滲濾液收集池池底的防滲層破損，滲濾液中的污染物砷以0.224 mg/L的濃度通過破損區(qū)域泄漏到土壤中，泄漏時(shí)長(zhǎng)是9125天，包氣帶厚度是9.0 m，土壤水中砷的初始濃度為0 mg/L，用HYDRUS-1D軟件建立數(shù)學(xué)模型，模擬砷在土壤水中的運(yùn)移規(guī)律。通過模擬結(jié)果可知：在僅考慮土壤水的對(duì)流和彌散的作用下，泄漏事故發(fā)生后，砷向下遷移形成垂向污染暈。隨時(shí)間增加，污染暈下滲深度逐漸增加，土壤水中砷的濃度逐漸增加。泄漏發(fā)生后第12.4天，潛水面處土壤水中砷濃度由0 mg/L增加至0.00001 mg/L，砷穿透整個(gè)包氣帶土層到達(dá)地下水中；泄漏后第9125天，土壤水中砷濃度達(dá)到峰值，為0.224 mg/L；此后，滲濾液停止泄漏，土壤水中砷濃度逐漸降低，至第10 220天，土壤水中砷濃度最大值為0.0362 mg/L。

通過模擬預(yù)測(cè)結(jié)果可知：土壤對(duì)污染物雖然具有一定的吸附阻截能力，但不能無限制地對(duì)污染物進(jìn)行阻截，在污染源持續(xù)泄漏的情況下，污染物會(huì)穿透土壤層進(jìn)入地下水中，進(jìn)而污染地下水。因此，地下水污染防治工作中還要重視土壤污染的防治。

該模型只考慮土壤中的對(duì)流、彌散作用對(duì)污染物擴(kuò)散的影響，未考慮土壤吸附、生物降解和化學(xué)反應(yīng)等作用對(duì)污染物轉(zhuǎn)化的影響，模擬預(yù)測(cè)的濃度值比實(shí)際值偏高，后續(xù)研究時(shí)應(yīng)盡量把上述各影響因素考慮在內(nèi)。另外，HYDRUS-1D軟件是針對(duì)污染物在土壤水中運(yùn)移規(guī)律的預(yù)測(cè)，預(yù)測(cè)結(jié)果表示污染物在土壤水中的濃度(單位是mg/L)，而非土壤自身的質(zhì)量濃度。但是《土壤環(huán)境質(zhì)量 建設(shè)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)(試行)》(GB 36600—2018)和《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)(試行)》(GB 15618—2018)中給出的管控標(biāo)準(zhǔn)值均為土壤的質(zhì)量濃度。因此，如果要對(duì)照上述標(biāo)準(zhǔn)對(duì)土壤質(zhì)量進(jìn)行評(píng)價(jià)，還需根據(jù)土壤容重和含水量等土壤理化性質(zhì)參數(shù)把預(yù)測(cè)結(jié)果換算成質(zhì)量濃度。