一維水分運(yùn)移模型水力參數(shù)敏感性分析

劉明明，段 磊，張 琛，宋 浩，呂婷婷，霍世璐

(1.長安大學(xué)環(huán)境科學(xué)與工程學(xué)院，西安 710054；2.旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室，西安 710054)

近年來，隨著計(jì)算的發(fā)展，利用數(shù)值模擬手段，研究包氣帶水分運(yùn)移已經(jīng)日益引起了人們的重視，并逐漸成為研究包氣帶水分運(yùn)移規(guī)律的重要手段[1-3]，但是，由于模型參數(shù)的不確定性影響了模型精確度，給建立包氣帶一維水分運(yùn)移模型的取值，造成了很大的困難,所以為了提高模型精度，針對參數(shù)的不確定性，進(jìn)行定量的分析與研究有著重大意義。

在數(shù)值模擬中，分析和測定模型參數(shù)對模型輸出結(jié)果的影響，這就是參數(shù)的敏感性分析[4,5]。目前，針對參數(shù)的不確定性對包氣帶水分運(yùn)移模型的研究比較多。王志濤等[6]利用HYDRUS-1D，采用單因素?cái)_動(dòng)法，得出n、θs和k的擾動(dòng)對垂直濕潤距離計(jì)算值和累積入滲量的影響最大；畢經(jīng)緯[7]、孟立江[8]以底部通量作為敏感性分析的指標(biāo)，得出飽和水力傳導(dǎo)度是最敏感參數(shù)；張?zhí)淼萚9]利用路面模型，得出淺層土壤水分對土壤質(zhì)地較為敏感，敏感性系數(shù)達(dá)到了0.45以上 ，并且砂土含量對土壤水分的影響更為顯著 ；但是，針對參數(shù)不確定性對非均質(zhì)剖面含水率空間分布的研究比較少。因此，本文基于HYDRUS[10]模型利用局部敏感性分析的方法，對模型參數(shù)進(jìn)行敏感性分析，旨在找出影響剖面含水率空間分布的顯著參數(shù)，給借助模型研究土壤水分動(dòng)態(tài)提供參考。

1 實(shí)驗(yàn)區(qū)概況

實(shí)驗(yàn)區(qū)位于新疆霍城縣伊犁河流域中下游的細(xì)土平原區(qū)，該地區(qū)屬溫帶半干旱氣候，平均氣溫8.2～9.4 ℃，年平均蒸發(fā)量約850～1 000 mm，年平均降水量約140～460 mm。

為觀測灌水之后剖面含水率的變化，利用含水率觀測儀器TDR，分別對剖面在0、10、20、40、60、90、140、190、200和240 cm深度的含水率進(jìn)行監(jiān)測，監(jiān)測頻率為次/3 h。

2 數(shù)值模擬

2.1 概念模型

根據(jù)實(shí)驗(yàn)場概況，建立一維水分運(yùn)移模型，模擬裸土試驗(yàn)場在灌水之后，水分在剖面上的空間分布。層狀非均勻包氣帶巖性如表1所示，其中主要巖性為粉質(zhì)黏土，并含有中砂以及粉土夾層(圖1)。因此，基于包氣帶巖層的考慮，本文將模型概化為四層。上邊界為開放邊界，接受降水、灌溉水補(bǔ)給及土壤蒸發(fā)排泄，蒸發(fā)量根據(jù)中國氣象網(wǎng)提供的氣象數(shù)據(jù)，利用Peman-Monteith計(jì)算得到[11,14]。潛水埋深為4.5 m，大于潛水面毛細(xì)上升高度，因此，本文將下邊界設(shè)置為自由排邊水邊界。

表1 土樣顆粒分析實(shí)驗(yàn)表Tab.1 Soil sample particle analysis experiment table

圖1 巖性剖面圖Fig.1 Lithology profile

2.2 數(shù)學(xué)模型

(1)

式中：h為土壤負(fù)壓，cm；c(h)為容水度，1/cm；h0為初始土壤負(fù)壓值，cm；P為降雨量，cm/d；I為灌溉量，cm/d；E為蒸發(fā)量，cm/d；t為時(shí)間，d；B為模型下邊界；

HYDRUS采用van-Genuchten模型描述土壤水分特征曲線和土壤導(dǎo)水率曲線兩類參數(shù)[15]，模型表達(dá)公式如下：

(2)

(3)

(4)

式中：θ(h)為土壤的特征曲線；θ為土壤含水率，%；θr為土壤的殘余含水率，%；θs為土壤的飽和含水率，%；h為壓力水頭；K(h)為壓力水頭下的導(dǎo)水率，%；ks為飽和導(dǎo)水率，%；α、n和m為van-Genuchten模型的形狀參數(shù)；Se為有效飽和度。

2.3 模型的離散化

采用等間隔的方法，以1 cm為間隔將包氣帶剖分為241個(gè)單元格，模擬期從2016年10月1日-2017年9月31日，從中分析參數(shù)對剖面含水率空間分布的影響，模型采用變時(shí)間間隔的方法，初始時(shí)間間隔為0.001 d，最小時(shí)間間隔為0.001 d，最大時(shí)間間隔為2 d。

2.4 土壤參數(shù)的識別

利用篩選法將用環(huán)刀在實(shí)驗(yàn)場取的原狀土樣進(jìn)行定名，再利用TST-70飽和滲透儀，測出飽和導(dǎo)水率，最后利用德國生產(chǎn)的ku-pF非飽和導(dǎo)水率測定系統(tǒng)[16]，得到負(fù)壓和含水率的關(guān)系，再利用MATLAB最小二乘法對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行參數(shù)擬合[17-19]，得到土壤水力參數(shù)如表2。

基于殘差平方和均方根誤差進(jìn)行誤差分析，并從表3的結(jié)果可知，土壤水分特征曲線擬合比較好，符合要求，因此本文將表2作為模型的初始值。

表2 土壤水分特征曲線擬合參數(shù)以及實(shí)測參數(shù)Tab.2 Soil moisture characteristic curve fitting parameters and measured parameters

表3 土壤水分特征曲線擬合誤差表Tab.3 Soil moisture characteristic curve fitting error table

2.5 模型識別與驗(yàn)證

模擬期內(nèi)經(jīng)過兩次灌水試驗(yàn)，因此本文分別將兩次灌水后的第8 d剖面含水率作為識別與驗(yàn)證期，在經(jīng)過對比研究后，發(fā)現(xiàn)擬合效果良好(表4、圖2)。

表4 識別后土壤水力參數(shù)Tab.4 Soil hydraulic parameters after identification

圖2 土壤剖面含水率擬合圖Fig.2

3 參數(shù)敏感性分析

3.1 參數(shù)敏感性分析理論依據(jù)

局部分析的方法，具有高效以及操作簡單的優(yōu)點(diǎn)。因此，本文將利用局部分析的方法對參數(shù)空間中的一點(diǎn)進(jìn)行敏感性分析[20]。

本文使用相對敏感性值將參數(shù)進(jìn)行歸一化處理，計(jì)算出敏感指數(shù) ，并根據(jù)敏感數(shù)進(jìn)行分類[21,22]。

I=(ΔO/ΔFi) (Fi/O)

(5)

式中：O為模型輸出結(jié)果；Fi為影響O的因子(參數(shù))；ΔO為模型輸出結(jié)果的改變值；ΔFi為影響O的因子(參數(shù))的改變量。敏感性分類見表5。

表5 敏感性分類Tab.5 Sensitivity classification

3.2 參數(shù)取值

根據(jù)上述巖層殘余含水率與室外實(shí)測初始含水率的對比發(fā)現(xiàn)，當(dāng)殘余含水率增加的值超過20%時(shí)，殘余含水率將高于巖層最低初始含水率，因此為了便于比較參數(shù)對剖面含水率的敏感性，本文將水力參數(shù)上下變幅統(tǒng)一變?yōu)?0%，以此來計(jì)算水力參數(shù)的敏感指數(shù)，表6為水力參數(shù)取值范圍。

表6 水力參數(shù)的取值Tab.6 Value of hydraulic parameters

3.3 土壤含水量模型參數(shù)的敏感性分析

本文利用參數(shù)局部敏感性分析的手段，通過改變水力參數(shù)的單一變量，來研究灌水15 cm之后，水力參數(shù)對剖面含水率影響的敏感指數(shù)，結(jié)果如表7所示。

表7 裸地情況下影響土壤含水量模型水力參數(shù)的敏感性分析Tab.7 Sensitivity analysis of hydraulic parameters affecting soil water content model under bare soil conditions

本文認(rèn)為當(dāng)敏感指數(shù)絕對值大于0.05時(shí)，即可認(rèn)為輸出結(jié)果對參數(shù)變化敏感(分級按照水力參數(shù)對剖面含水率影響的最大敏感指數(shù)進(jìn)行歸類)。

從表7參數(shù)的敏感指數(shù)可知，粉質(zhì)黏土層θs1上浮20%時(shí)，敏感指數(shù)屬于Ⅲ級，最大敏感指數(shù)最大達(dá)到0.60，對土壤上部粉質(zhì)黏土剖面含水率，影響比較大，但是對下部其他巖層的影響比較小，當(dāng)θs1下浮20%時(shí)，敏感指數(shù)屬于Ⅲ級，最大敏感指數(shù)達(dá)到0.91，但影響范圍只局限于本巖層。中砂巖層θr2上下浮動(dòng)20%時(shí)，敏感指數(shù)幾乎相等，同時(shí)達(dá)到Ⅲ級，最大敏感指數(shù)達(dá)到0.98，對本層的剖面含水率影響最大，對其他巖層含水率并沒有造成太大的影響。粉土層 上浮和下浮20%時(shí)，敏感指數(shù)差不多，最大達(dá)到0.72，同屬Ⅲ級，對第三層剖面含水率影響比較大，當(dāng)n3上下浮動(dòng)時(shí)，發(fā)現(xiàn)當(dāng)n3向下浮動(dòng)時(shí)，其敏感指數(shù)達(dá)到了Ⅳ級，說明了n3是對第三層剖面含水率影響最大的參數(shù)。當(dāng)θs4上浮和下浮20%時(shí)，敏感指數(shù)都達(dá)到Ⅲ級，但n4對剖面含水率影響指數(shù)達(dá)到了Ⅳ級，表明了在下部粉質(zhì)黏土層中，n4對該層剖面含水率影響最大，而且對上部粉土層也產(chǎn)生了很大的影響。

為了驗(yàn)證上述敏感性分析結(jié)果，本次模擬將灌水變?yōu)?0 cm，來研究水力參數(shù)對剖面含水率的影響(表8)。

表8 裸地情況下影響土壤含水量模型水力參數(shù)的敏感性分析Tab.8 Sensitivity analysis of hydraulic parameters affecting soil water content model under bare soil conditions

從表8參數(shù)的敏感指數(shù)可知，上部粉質(zhì)黏土n1和θs1對剖面含水率的影響最大，不過n1對該層剖面含水率最大，而對于中砂層來說，θr2對該層的剖面含水率，起著決定性作用，而其他水力參數(shù)，對剖面含水率幾乎沒有造成影響，而粉土和下部粉質(zhì)黏土層的n和θs共同影響著剖面含水率，但是n對該層剖面含水率最大，由此不難看出，灌水15 cm和灌水20 cm時(shí)，得出的結(jié)論一致。

3.4 分析與討論

通過對表7和表8的對比發(fā)現(xiàn)，在上部粉質(zhì)黏土層中，n1和θs1對第一層剖面含水率敏感指數(shù)最高，并且敏感指數(shù)隨著深度的增加而減少。在中砂層中，只有θr2對本層剖面含水率敏感性最大，其他水力參數(shù)上下浮動(dòng)，難以對剖面含水率產(chǎn)生很大的影響，粉土巖層以及下部粉質(zhì)黏土層中n和θs對本層的含水率，有著重大的影響，由此可以看出，在細(xì)顆粒中n和θs對剖面含水率影響比較大，而且n的敏感指數(shù)要大于θs的敏感指數(shù)，這主要是VG模型參數(shù)中θs對水分特征曲線低吸引力段影響顯著，而n變化會影響土壤水分特征曲線的彎曲程度以及整體形狀，但其參數(shù)主要在土壤釋水的后期，發(fā)揮著主導(dǎo)性作用[23,24]，因而在灌水后期，細(xì)顆粒巖層n對剖面含水率的影響程度大于θs。VG模型中θs在灌水后期開始發(fā)揮著重要作用[25]，而中砂層由于顆粒比較粗，所以土壤水在前期釋水的過程中，θs2和n2發(fā)揮著重要作用，但是隨著釋水過程的推進(jìn)，延遲與滯緩的時(shí)間，對剖面含水率影響逐漸減弱，當(dāng)剖面含水率降低接近該巖層θr2時(shí)，θr2對該剖面含水率影響占據(jù)了主導(dǎo)作用。

4 結(jié) 語

水分經(jīng)過300多天的水分運(yùn)移之后， 在非均質(zhì)巖層細(xì)顆粒以及粗顆粒巖層中，α和k參數(shù)的變化難以對剖面含水率產(chǎn)生太大的影響，而在細(xì)顆粒巖層中，θs和n對剖面含水率發(fā)揮著重要著作用，在粗顆粒巖層中，θr對剖面含水率的影響最大。

數(shù)值模型在實(shí)際應(yīng)用的過程中，由于巖層復(fù)雜多變，給模型參數(shù)的取值造成了很大困難，因此本次研究結(jié)論，可以減少在非均質(zhì)中構(gòu)建模型時(shí)，需要調(diào)整的參數(shù)，對于提高模型精確度有著重大意義，所以，在進(jìn)行數(shù)值模擬時(shí)，應(yīng)對細(xì)顆粒巖層的θs進(jìn)行加密測量，而對n要進(jìn)行重點(diǎn)調(diào)參，在粗顆粒的巖層中，要提高對θr的重視程度。