高志鵬，屈吉鴻，陳南祥，李五金(華北水利水電大學(xué)資源與環(huán)境學(xué)院，鄭州 450045)

在自然界中，溶質(zhì)運移在飽和、非飽和含水系統(tǒng)中廣泛存在，非飽和帶的水力參數(shù)對水分和溶質(zhì)的運移具有較大的影響[1]。因此，研究非飽和土壤的水力參數(shù)對非飽和、非飽和-飽和含水系統(tǒng)中污染物的入滲影響是至關(guān)重要的，學(xué)者們通過大量理論和實驗研究發(fā)現(xiàn)了土壤水分特征曲線，發(fā)現(xiàn)土壤水勢與含水率之間存在著一定的規(guī)律性，但是不管是室內(nèi)實驗還是野外田間試驗，測定土壤水分特征曲線均存在一定的局限性，為解決這一問題，國內(nèi)外學(xué)者為此建立了經(jīng)驗與半經(jīng)驗理論模型[2]，其中應(yīng)用最廣的是van Genuchten提出的計算土壤含水率的經(jīng)驗公式。非飽和土壤的水力參數(shù)控制含水率的分布，進而影響非飽和帶溶質(zhì)運移的過程[3]。

如今，采用VG模型研究非飽和帶的水分和溶質(zhì)運移的研究較多[4-8]，而非飽和土壤的水力參數(shù)對溶質(zhì)在非飽和含水系統(tǒng)中遷移規(guī)律影響的研究卻比較少。Liu等通過研究變密度流下溶質(zhì)在非飽和-飽和含水系統(tǒng)中遷移的規(guī)律，探討了非飽和帶主要水力參數(shù)(α、n、θr)對溶質(zhì)穿透非飽和-飽和過渡帶的影響，指出低含水率、中含水率、高含水率具有不同的影響規(guī)律[1]。Pan等采用基于方差分解的回歸分析方法(The sample-based regression and decomposition methods)對影響非飽和帶中水流及溶質(zhì)運移的參數(shù)進行敏感性及相關(guān)性分析，研究分別對參數(shù)進行局部敏感性分析及全局敏感性分析，通過判定各個輸入?yún)?shù)(van Genuchtenα、n等)對輸出變量(滲流量、溶質(zhì)通量等)的方差貢獻率的大小，進而估計各個輸入?yún)?shù)的重要程度[9]。王志濤等通過構(gòu)建一維垂直入滲模型，采用單因素擾動分析法研究了VG模型中5個參數(shù)的擾動對垂直濕潤距離和累計入滲量的影響[10]。

對非飽和帶水分及溶質(zhì)運移開展參數(shù)敏感性分析對于水資源管理、環(huán)境保護以及修復(fù)模式的提出均具有重大的意義。敏感性分析研究的是輸入?yún)?shù)對輸出結(jié)果的影響[11]，包括局部敏感性分析和全局敏感性分析，局部敏感性分析評估了輸出結(jié)果隨某一個參數(shù)變化的過程，而全局敏感性分析則評估了輸出結(jié)果隨著多種輸入?yún)?shù)變化而改變的過程，局部敏感性分析能直觀的展現(xiàn)各參數(shù)對結(jié)果的影響程度，區(qū)分出重要參數(shù)及相對不重要參數(shù)，但是忽略了參數(shù)之間的相互影響；而全局敏感性分析則在考慮各參數(shù)間相互影響的條件下研究對輸出結(jié)果的影響。目前國內(nèi)對地下水?dāng)?shù)值模型的參數(shù)敏感性分析大都采用局部敏感性分析，典型做法是只改變一個參數(shù)值，觀察其對模擬結(jié)果的影響程度，進而計算敏感度[12]。

衛(wèi)河流域地處海河流域的南部，河南省的北部。衛(wèi)河河水及沿岸地下水是流域內(nèi)重要的淡水資源，水質(zhì)的變化將影響到流域內(nèi)人類生活和生態(tài)環(huán)境的穩(wěn)定。衛(wèi)河常年接收工業(yè)廢水和生活污水，嚴(yán)重污染的河水通過包氣帶自由滲漏勢必影響到地下水的水質(zhì)。包氣帶是河水補給地下水必經(jīng)的通道，因此研究包氣帶水力參數(shù)的擾動對水流及溶質(zhì)運移的影響可以為衛(wèi)河流域地下水保護、修復(fù)和治理提供理論基礎(chǔ)。

本文以描述非飽和土壤水力特征的VG模型為理論基礎(chǔ)，采用Hydrus-1d構(gòu)建一維非飽和流數(shù)值模型，通過構(gòu)建不同的情景模式，研究壓力水頭和NaCl濃度隨深度的變化規(guī)律，探討了VG模型參數(shù)對水流及溶質(zhì)運移變化的敏感度大小。

1 模型的建立

1.1 Hydrus-1d模型簡介

Hydrus-1d是美國農(nóng)業(yè)部、美國鹽堿實驗室等機構(gòu)開發(fā)出來的[13]，是一種基于有限元的，可用來模擬一維變飽和度地下水流、根系吸水、溶質(zhì)運移和熱運移的計算機模型。模型可以處理各類復(fù)雜邊界，包括定水頭和變水頭邊界、定水頭梯度邊界、定流量邊界、滲水邊界、自由排水邊界、大氣邊界等[14]，同時具有靈活的輸入輸出功能，在土壤中水分運動、鹽分、農(nóng)藥及土壤氮素運移方面得到了廣泛應(yīng)用[15]。

1.2 模型設(shè)置

模型模擬地下0～100 cm深度范圍內(nèi)的土壤，剖分1層，土壤類型采用軟件提供的壤土(Loam)，模擬時段為1 d，采用變時間步長輸出模擬結(jié)果(0.4、0.8、1 d)。土壤水流模型采用van Genuchten-Mualem，不考慮水分滯后效應(yīng)。水流模型上邊界為定水頭邊界，邊界水頭值取為1 cm，下邊界為自由排水邊界。入滲過程中，濕潤鋒未達(dá)到下邊界[20]，故下邊界可概化為定水頭或隔水邊界。溶質(zhì)運移上邊界為定濃度邊界，為10 mg/L，下邊界為零濃度梯度邊界。

1.3 土壤水力函數(shù)方程

Hydrus-1d的土壤水力函數(shù)方程為van Genuchten-Mualem公式[16]，其表達(dá)形式為：

K(h)=krKs

(1)

kr=θ1/2e[1-(1-θ1/me)m]2

(2)

(3)

θ=θr+θe(θs-θr)

(4)

(5)

式中：h為壓力水頭；kr為相對滲透系數(shù)；Ks為飽和導(dǎo)水率；θ為含水率；θs為飽和含水率；θr為殘余含水率；θe為有效含水率；α、n為經(jīng)驗擬合參數(shù)；h<0表示包氣帶，h≥0表示飽和帶。

1.4 水流模型

Hydrus-1d軟件模擬土壤水分運動的方程為Richards方程[17]，構(gòu)建一維垂向均質(zhì)非飽和入滲水流模型為：

(6)

式中：C(h)=dθ/dh，為容水度；K(h)為非飽和導(dǎo)水率；h為壓力水頭；z為垂向坐標(biāo)，規(guī)定z向上為正；t為入滲時間；h0(z)為初始壓力水頭分布；L為入滲深度。

1.5 溶質(zhì)運移模型

Cl-是公認(rèn)的保守性離子，其在地下水中具有高溶解性而又難以被植物、細(xì)菌或土壤介質(zhì)所吸收[18]，因此選擇NaCl為溶質(zhì)運移的模擬對象。采用傳統(tǒng)的對流-彌散方程(Convection-Dispersion Equation, CDE)來描述溶質(zhì)運移的過程[19]。一維垂向溶質(zhì)運移模型可描述為：

(7)

式中：C為NaCl濃度；θ為含水率；DL為縱向彌散系數(shù)；qz為垂向水分達(dá)西流速；L為運移深度。

1.6 情景設(shè)置

由式(1)～(5)可知，VG模型中共有θr、θs、α、n、Ks5個參數(shù)，以軟件提供的壤土的經(jīng)驗參數(shù)值為基準(zhǔn)參數(shù)(STD)構(gòu)建模型，采用擾動分析法[21]，將各參數(shù)在基準(zhǔn)參數(shù)的基礎(chǔ)上分別上下擾動5%、10%。模擬時，將某一個參數(shù)作一定的擾動而不改變其他參數(shù)，重新運行模型，得到各情景下壓力水頭和溶質(zhì)濃度隨深度的變化值。各情景設(shè)置見表1。

表1 模擬情景設(shè)置Tab.1 Situation modes

2 敏感性分析

2.1 敏感系數(shù)

敏感度分析是衡量改變一個因子時對其他變量影響的量，進行敏感度分析是為了確定模型的結(jié)果對模型參數(shù)的敏感程度，它將提供有關(guān)模型參數(shù)如何影響模擬結(jié)果不確定性的重要信息。如果模擬結(jié)果對某一特定參數(shù)高度敏感，那么模型做出重要解釋和預(yù)報的能力將受到和該參數(shù)有關(guān)的不確定性的嚴(yán)重影響[22]。模型因變量對一個輸入?yún)?shù)的敏感性是該因變量對該參數(shù)的偏微分[23]：

(8)

方程進一步標(biāo)準(zhǔn)化為量綱為一的形式為：

(9)

敏感度分析時，研究參數(shù)(如飽和滲透系數(shù)Ks)的變化多大(即飽和滲透系數(shù)的改變ΔKs/K0s)會引起數(shù)值解多大的變化(如模型變量壓力水頭的變化Δh)。給出某參數(shù)不同的擾動值(如ΔKs/K0s)，通過計算壓力水頭的變化率Δh/h0或溶質(zhì)濃度的變化率Δc/c0，然后繪制ΔKs/K0s和Δh/h0(Δc/c0)的關(guān)系曲線，則曲線的斜率即為敏感系數(shù)Xi,k。式(9)中Xi,k為正值時，表示y隨ak的增大而增大；Xi,k為負(fù)時，表示表示y隨ak的增大而減小。Xi,k大于1(或小于-1)，表明輸出變量的變化較研究參數(shù)的變化大，說明該參數(shù)對輸出變量的影響較大。

2.2 敏感度指數(shù)

Xi,k的絕對值表示模型輸入?yún)?shù)對模型輸出變量的相對敏感性，或稱敏感度指數(shù)。定義敏感度指數(shù)L為：

(11)

式中：m為觀測點總數(shù)，本次取10；n為擾動次數(shù)，本次取4；Δi為擾動值(分別為-10%、-5%、5%、10%)。

上式定義的敏感度指數(shù)L為平均值，可以定量求得參數(shù)的敏感性，但是具有一定的粗糙性。

3 結(jié)果與分析

3.1 VG模型參數(shù)的擾動對水流及溶質(zhì)運移的影響

3.1.1 θr擾動的影響

圖1(a)和圖1(b)為θr的擾動對壓力水頭的影響，圖1分別表示了3個時刻(0.4 d、0.8 d、1.0 d)各擾動情景下壓力水頭隨深度的變化規(guī)律。由圖1(a)可以看出，θr的擾動對壓力水頭的影響較小(觀察同一模擬期壓力水頭曲線所圍面積，面積越大，影響越大)，對濕潤鋒入滲深度的影響也很小，幾乎可以忽略[22]；從圖1(b)可以看出，θr對壓力水頭的影響規(guī)律接近線性正相關(guān)，即敏感系數(shù)為正值，說明θr增大，壓力水頭也增大，并且隨著入滲時間的推移，其影響程度呈先增大后減小的趨勢；斜率小于1說明θr的影響較小。

圖1(c)和圖1(d)為θr的擾動對溶質(zhì)運移的影響。從圖1(c)可以看出，溶質(zhì)濃度隨深度的變化曲線在A1～A4情景下相對基準(zhǔn)情景(STD)基本無變化，說明θr對溶質(zhì)運移的影響非常小，幾乎可以忽略；從圖1(d)可以看出，其影響規(guī)律也接近線性正相關(guān)。

圖1 θr擾動的影響Fig.1 Influence of disturbance of θr

3.1.2 θs擾動的影響

圖2(a)和圖2(b)為θs的擾動對壓力水頭的影響。由圖2(a)可以看出θs對壓力水頭的影響與θr相比更大；B3、B4情景下θs較基準(zhǔn)情景減小(即負(fù)擾動)，濕潤鋒較基準(zhǔn)情景(STD)下滲深度大，θs增大，B1、B2濕潤鋒較基準(zhǔn)情景(STD)下滲深度小，說明θs對濕潤鋒具有較大的影響[10]；其影響規(guī)律為非線性負(fù)相關(guān)，斜率為負(fù)，說明壓力水頭或濕潤鋒入滲深度隨θs的增大而減?。磺译S著入滲時間的推移，其影響程度逐漸增大[圖2(b)]。通常，飽和含水率(θs)與土壤的孔隙度相等，而土壤孔隙的多少，決定了巖土儲水容水的能力，而孔隙的大小，對于巖土滯留、釋出及傳輸水的能力影響很大[24]，改變土壤的飽和含水率，相當(dāng)于改變了土壤的儲水、導(dǎo)水的能力，因此，對土壤的飽和含水率(θs)施加一定的擾動，對濕潤鋒的前進以及壓力水頭的影響會很大。

圖2(c)和圖2(d)為θs對溶質(zhì)運移的影響。隨著入滲時間的推移，土壤層上部飽和帶的厚度逐漸增加(如0.4 d時0～20 cm深度，0.8 d時0～40 cm深度為飽和帶)，在濕潤鋒前進的前部則形成非飽和帶，圖2(c)表明θs的擾動對飽和帶與非飽和帶溶質(zhì)的濃度隨深度的變化均具有較大影響。從圖2(d)可以看出θs對溶質(zhì)運移的影響呈非線性負(fù)相關(guān)，且隨著時間的推移，其影響呈先增大后減小的趨勢。與對壓力水頭的影響相比，θs對溶質(zhì)運移的影響正好相反，其負(fù)擾動的影響很大，而正擾動的影響則相對較小。

圖2 θs擾動的影響Fig.2 Influence of disturbance of θs

3.1.3 α擾動的影響

圖3(a)和圖3(b)為α的擾動對壓力水頭的影響。從圖3(a)可以看出α對壓力水頭和濕潤鋒的影響也較大，α的負(fù)擾動越大，濕潤鋒的下滲深度越大，α的正擾動越大，濕潤鋒的前進深度比基準(zhǔn)情景的深度小10 cm左右，當(dāng)α的擾動為-10%時，1 d時濕潤鋒甚至運移到了模型最底部；α對壓力水頭的影響也比較大，且其影響規(guī)律呈非線性負(fù)相關(guān)，與θs相同，α對壓力水頭的敏感系數(shù)也為負(fù)，說明α的增大伴隨著壓力水頭和濕潤鋒入滲深度的減小[圖3(b)]。

圖3(c)和圖3(d)為α的擾動對溶質(zhì)運移的影響。由圖3(c)可以看出α對溶質(zhì)運移的影響相對較小，α的擾動不改變?nèi)苜|(zhì)濃度隨深度變化的整體趨勢，圖3(d)表明α對溶質(zhì)運移的影響規(guī)律為負(fù)相關(guān)，且負(fù)擾動為線性關(guān)系，正擾動為非線性關(guān)系。其敏感系數(shù)為負(fù)同樣說明α的增大會引起溶質(zhì)濃度的減小。

圖3 α擾動的影響Fig.3 Influence of disturbance of α

3.1.4 n擾動的影響

圖4(a)和圖4(b)為n的擾動對壓力水頭的影響。從圖4(a)可以看出，n對濕潤鋒和壓力水頭的影響也比較大，與α相反，n的正擾動會加快濕潤鋒的前進，負(fù)擾動則會減緩。圖4(b)說明n對壓力水頭的影響規(guī)律為非線性正相關(guān)，隨著入滲時間的推移其影響逐漸變大，負(fù)擾動對壓力水頭的影響比正擾動的影響大的多。

圖4(c)和圖4(d)為n的擾動對溶質(zhì)運移的影響。從圖4可以看出n對溶質(zhì)運移的影響很大，其影響為非線性正相關(guān)。其對溶質(zhì)濃度的影響與壓力水頭的影響正好相反，輕微的正擾動可能引起相當(dāng)大的溶質(zhì)濃度的波動，負(fù)擾動的影響則相對較小。

3.1.5 Ks擾動的影響

圖5(a)和圖5(b)為Ks的擾動對壓力水頭的影響。飽和滲透系數(shù)(Ks)可以定量說明土壤的滲透性能，滲透系數(shù)越大，巖石的滲透能力越強，圖5(a)中飽和滲透系數(shù)越大(如1d-E1<1d-E2<1d-STD<1d-E3<1d-E4)，相同深度處的壓力水頭越大，濕潤鋒下滲的深度也越大。Ks的擾動對壓力水頭的影響呈非線性正相關(guān)，且隨著時間的推移，其影響程度也越來越大[圖5(b)]。

圖5(c)和圖5(d)為Ks的擾動對溶質(zhì)運移的影響。飽和滲透系數(shù)不僅對壓力水頭有較大的影響，因為其控制對流-彌散方程中的對流項，因此對溶質(zhì)運移也具有較大的影響，如圖5(c)和圖5(d)所示，其影響規(guī)律為非線性正相關(guān)，與壓力水頭不同，Ks的正擾動的影響程度更大。

3.2 敏感性分析

分別計算0.4 d、0.8 d、1.0 d時各參數(shù)對壓力水頭和溶質(zhì)濃度的平均敏感度指數(shù)L，結(jié)果見圖6。從圖6可以看出各參數(shù)對壓力水頭的敏感度指數(shù)大小關(guān)系為：θr<α

圖4 n擾動的影響Fig.4 Influence of disturbance of n

圖5 Ks擾動的影響Fig.5 Influence of disturbance of Ks

圖6 壓力水頭和溶質(zhì)濃度對5個參數(shù)的敏感度Fig.6 The sensitivity of the pressure head and the solute concentration to the 5 parameters

4 結(jié) 語

VG模型中θr、θs、α、n、Ks5個參數(shù)對模型輸出變量的影響不盡相同，5個參數(shù)的擾動對壓力水頭和溶質(zhì)運移的影響不僅呈現(xiàn)出不同的規(guī)律，其影響程度也隨著模擬期的改變而變化。

(1)壓力水頭對5個參數(shù)的敏感度大小關(guān)系為：θr<α

(2)溶質(zhì)運移對5個參數(shù)的敏感度大小關(guān)系為：θr<α

(3)據(jù)以上分析，采用數(shù)值模型研究污染河水對溶質(zhì)運移的影響，應(yīng)保證模型中敏感度高的水力參數(shù)的精度。θs、Ks、n的敏感度高且隨模擬時間的增長對壓力水頭的影響呈上升趨勢，α和θr的敏感度相對較低且對壓力水頭和溶質(zhì)運移的影響均呈先上升后下降的趨勢，故應(yīng)保證θs、Ks、n的精度[10]。

[1] Liu Y，Kuang X X，Jiao J J，et al. Numerical study of variable-density flow and unsaturated-saturated porous media[J]. Journal of Contaminant Hydrology，2015，173：117-130.

[2] 錢天偉，劉春國. 飽和-非飽和土壤污染物運移[M]. 北京：中國環(huán)境科學(xué)出版社，2007.

[3] van Genuchten M T. A Closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of American Journal，1980，44(5)：892-898.

[4] Berlin M，Kumar G S，Nambi I M. Numerical modelling on transport of nitrogen from wastewater and fertilizer applied on paddy fields [J]. Ecological Modelling，2014，278：85-99.

[5] Berg S J，Illman W A. Improved predictions of saturated and unsaturated zone drawdowns in a heterogeneous unconfined aquifer via transient hydraulic tomography：laboratory sandbox experiments[J]. Journal of Hydrology，2012，470-471：172-183.

[6] 劉 超，何江濤，沈 楊，等. 非均質(zhì)包氣帶三氮累積轉(zhuǎn)化模擬砂箱試驗[J]. 地學(xué)前緣，2012,19(6)：236-242.

[7] 龐雅婕，劉長禮，王翠玲，等. 某化工廠廢液CODCr在包氣帶中的遷移規(guī)律及數(shù)值模擬[J]. 水文地質(zhì)工程地質(zhì)，2013,40(3)：115-120.

[8] 丁素玲. 非均質(zhì)包氣帶中“三氮”遷移轉(zhuǎn)化數(shù)值模擬[D]. 北京：中國地質(zhì)大學(xué)，2012.

[9] Pan F，Zhu J T，Ye M，et al. Sensitivity analysis of unsaturated flow and contaminant transport with correlated parameters[J]. Journal of Hydrology，2011，397：238-249.

[10] 王志濤，繳錫云，韓紅亮，等. 土壤垂直一維入滲對VG模型參數(shù)的敏感性分析[J]. 河海大學(xué)學(xué)報(自然科學(xué)版)，2013,41(1)：80-84.

[11] Jacques J, Lavergne C, Devictor N. Sensitivity analysis in presence of model uncertainty and correlated inputs[J]. Reliability Engineering & System Safety ， 2006，91：1 126-1 134.

[12] 李木子. 地下水環(huán)境影響評價數(shù)值模擬中的參數(shù)敏感性分析[D]. 北京：中國地質(zhì)大學(xué)，2013.

[13] 李 瑋，何江濤，劉麗雅，等. Hydrus-1d軟件在地下水污染風(fēng)險評價中的應(yīng)用[J]. 中國環(huán)境科學(xué)，2013,33(4)：639-647.

[14] 余根堅，黃介生，高占義. 基于Hydrus模型不同灌水模式下土壤水鹽運移模擬[J]. 水利學(xué)報，2013，44(7)：826-834.

[15] 郝芳華，孫 雯，曾阿妍，等. Hydrus-1d模型對河套灌區(qū)不同施灌情景下氮素遷移的模擬[J]. 環(huán)境科學(xué)學(xué)報，2008，28(5)：853-858.

[16] Hsin-Chi J Lin，David R Richards ，Cary A Talbot，et al. FEMWATER：a three-dimensional finite element computer model for simulating density-dependent flow and transport in Variably Saturated Media Version 3.0[Z].

[17] ZHU Y，ZHA Y Y，TONG J X，et al. Method of coupling 1-D unsaturated flow with 3-D saturated flow on large scale[J]. Water Science and Engineering，2011,4(4)：353-373.

[18] 沈照理，朱宛華，鐘佐燊. 水文地球化學(xué)基礎(chǔ)[M]. 北京：地質(zhì)出版社，1993.

[19] 王小丹，鳳 蔚，王文科，等. 基于 HYDRUS-1D模型模擬關(guān)中盆地氮在包氣帶中的遷移轉(zhuǎn)化規(guī)律[J]. 地質(zhì)調(diào)查與研究，2015，38(4)：291-298,304.

[20] 范嚴(yán)偉，趙文舉，畢貴權(quán). van Genuchten模型參數(shù)變化對土壤入滲特性的影響分析[J]. 中國農(nóng)村水利水電，2016,(3)：52-56.

[21] 魯程鵬，束龍倉，劉麗紅，等. 基于靈敏度分析的地下水?dāng)?shù)值模擬經(jīng)度適應(yīng)性評價[J]. 河海大學(xué)學(xué)報(自然科學(xué)版)，2010，38(1)：26-30.

[22] 薛禹群，謝春紅. 地下水?dāng)?shù)值模擬[M]. 北京：科學(xué)出版社，2007.

[23] Zheng C M，Bennett G D. 地下水污染物遷移模擬[M]. 北京：高等教育出版社，2009.

[24] 張人權(quán)，梁 杏，靳孟貴，等. 水文地質(zhì)學(xué)基礎(chǔ)[M]. 北京：地質(zhì)出版社，2011.