一維水流及溶質運移對VG模型參數的敏感性分析

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

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

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

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

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

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

1 模型的建立

1.1 Hydrus-1d模型簡介

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

1.2 模型設置

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

1.3 土壤水力函數方程

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

K(h)=krKs

(1)

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

(2)

(3)

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

(4)

(5)

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

1.4 水流模型

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

(6)

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

1.5 溶質運移模型

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

(7)

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

1.6 情景設置

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

表1 模擬情景設置Tab.1 Situation modes

2 敏感性分析

2.1 敏感系數

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

(8)

方程進一步標準化為量綱為一的形式為：

(9)

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

2.2 敏感度指數

Xi,k的絕對值表示模型輸入參數對模型輸出變量的相對敏感性，或稱敏感度指數。定義敏感度指數L為：

(11)

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

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

3 結果與分析

3.1 VG模型參數的擾動對水流及溶質運移的影響

3.1.1 θr擾動的影響

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

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

圖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較基準情景減小(即負擾動)，濕潤鋒較基準情景(STD)下滲深度大，θs增大，B1、B2濕潤鋒較基準情景(STD)下滲深度小，說明θs對濕潤鋒具有較大的影響[10]；其影響規律為非線性負相關，斜率為負，說明壓力水頭或濕潤鋒入滲深度隨θs的增大而減小；且隨著入滲時間的推移，其影響程度逐漸增大[圖2(b)]。通常，飽和含水率(θs)與土壤的孔隙度相等，而土壤孔隙的多少，決定了巖土儲水容水的能力，而孔隙的大小，對于巖土滯留、釋出及傳輸水的能力影響很大[24]，改變土壤的飽和含水率，相當于改變了土壤的儲水、導水的能力，因此，對土壤的飽和含水率(θs)施加一定的擾動，對濕潤鋒的前進以及壓力水頭的影響會很大。

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

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

3.1.3 α擾動的影響

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

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

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

3.1.4 n擾動的影響

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

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

3.1.5 Ks擾動的影響

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

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

3.2 敏感性分析

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

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

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

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

4 結 語

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

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

(2)溶質運移對5個參數的敏感度大小關系為：θr<α

(3)據以上分析，采用數值模型研究污染河水對溶質運移的影響，應保證模型中敏感度高的水力參數的精度。θs、Ks、n的敏感度高且隨模擬時間的增長對壓力水頭的影響呈上升趨勢，α和θr的敏感度相對較低且對壓力水頭和溶質運移的影響均呈先上升后下降的趨勢，故應保證θ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]. 北京：中國環境科學出版社，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] 劉 超，何江濤，沈 楊，等. 非均質包氣帶三氮累積轉化模擬砂箱試驗[J]. 地學前緣，2012,19(6)：236-242.

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

[8] 丁素玲. 非均質包氣帶中“三氮”遷移轉化數值模擬[D]. 北京：中國地質大學，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模型參數的敏感性分析[J]. 河海大學學報(自然科學版)，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] 李木子. 地下水環境影響評價數值模擬中的參數敏感性分析[D]. 北京：中國地質大學，2013.

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

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

[15] 郝芳華，孫 雯，曾阿妍，等. Hydrus-1d模型對河套灌區不同施灌情景下氮素遷移的模擬[J]. 環境科學學報，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] 沈照理，朱宛華，鐘佐燊. 水文地球化學基礎[M]. 北京：地質出版社，1993.

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

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

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

[22] 薛禹群，謝春紅. 地下水數值模擬[M]. 北京：科學出版社，2007.

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

[24] 張人權，梁 杏，靳孟貴，等. 水文地質學基礎[M]. 北京：地質出版社，2011.