考慮蒸發(fā)鋒面的包氣帶一維液相運動方程解析解

劉夢茹, 程大偉, 盧玉東

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

1 研究背景

包氣帶是大氣水、地表水、地下水三者之間發(fā)生水力聯(lián)系并進行水分交換的場所，同時也是地表污染物進入地下的通道。因此包氣帶中液相運動的研究在地下水污染管理、維持地下水采補平衡[1]等領(lǐng)域中有著廣泛的應(yīng)用。包氣帶中土壤濕度是控制水分交換的關(guān)鍵因素，其狀態(tài)受季節(jié)和氣候變化的影響[2]。在雨季，降雨作用在包氣帶中形成的液相運動為穩(wěn)態(tài)入滲[3]；在旱季，土面蒸發(fā)和植物蒸騰是控制包氣帶中潛水位以上區(qū)域液相運動的主要因素。多數(shù)學者在推導(dǎo)包氣帶液相運動方程時僅考慮入滲情況，而由蒸發(fā)所引起的包氣帶液相運動在實際中大致可以分為以下兩種形式[4-7]：當?shù)叵滤簧疃茸銐驕\時，包氣帶液相運動完全受毛細作用控制，始終維持水力連續(xù)，且服從達西定律;當?shù)叵滤簧疃瘸^某一臨界深度時，受到重力因素及黏滯力的影響，蒸發(fā)鋒面[7]將會由地表向地下水位發(fā)展。

在推導(dǎo)包氣帶一維液相運動方程解析解過程中常用冪函數(shù)型或指數(shù)型滲透函數(shù)，諸如線性模型[8]、基于Brooks-Corey土水特征曲線建立的滲透性函數(shù)模型[9]、低飽和度下Van Genuchten 模型[10-11]、常數(shù)滲透系數(shù)模型[12]等均為冪函數(shù)型滲透性函數(shù)[3]。由于冪函數(shù)型滲透性函數(shù)在推導(dǎo)液相運動方程解析解時變量難以合并,導(dǎo)致所獲得的解析解常為負壓水頭的隱式表達，此外，在求解的過程中這些解僅考慮入滲情況。Gardner模型為指數(shù)型滲透系數(shù)函數(shù)，可對液相運動方程進行線性化處理，易推導(dǎo)解析解。在穩(wěn)態(tài)入滲或蒸發(fā)條件下,一維液相運動方程的解析解相比數(shù)值解簡潔[13],并且易基于相應(yīng)解析解建立反演模型,進行參數(shù)反演，從而具有廣泛的實際應(yīng)用,如通過實測數(shù)據(jù)反演估計水力參數(shù)、為更復(fù)雜的液相運動提供標準測試方案等。

本文基于指數(shù)型滲透系數(shù)模型Gardner模型求解擴散度，分別推導(dǎo)在穩(wěn)態(tài)入滲和穩(wěn)態(tài)蒸發(fā)兩種情況下包氣帶中均質(zhì)土和成層土液相運動方程的解析表達以及蒸發(fā)鋒面深度表達式。分析各水力參數(shù)對一維液相運動方程的影響并通過解析解進行水力參數(shù)的反演。

2 均質(zhì)土水力連續(xù)區(qū)有效飽和度剖面及蒸發(fā)鋒面深度的解答

LU Ning等[14]通過給定適當?shù)某跏紬l件和邊界條件，推導(dǎo)出以基質(zhì)吸力為變量的液相運動方程解析解。本節(jié)在此基礎(chǔ)上，推導(dǎo)出有效飽和度剖面解析解。

2.1 均質(zhì)土水力連續(xù)區(qū)液相運動方程解析解

一維穩(wěn)態(tài)流下包氣帶液相運動方程邊值問題可寫作如下形式：

(1)

式中：θ為含水量，m3/m3；D(θ)為擴散率；K(θ)為非飽和土滲透系數(shù)，m/s。

非飽和土的滲透系數(shù)是負壓水頭的函數(shù)，常被定義為：

K=KsKr

(2)

式中：Ks為飽和土滲透系數(shù)；Kr為相對滲透系數(shù)。

在非飽和土中,相對滲透系數(shù)Kr和含水量θ均為負壓水頭的函數(shù)[15]。Gardner模型常被用于求解析解，所描述的相對滲透系數(shù)與負壓水頭之間呈指數(shù)關(guān)系，即

Kr=eαgh

(3)

式中：αg為Gardner模型參數(shù)。

(4)

式中:Se為有效飽和度;θr為殘余含水量;θs為飽和含水量。公式(4)即為Srivastava提出的指數(shù)型土水特征關(guān)系[16]。

將公式(4)轉(zhuǎn)化為負壓水頭與有效飽和度的關(guān)系，則有：

(5)

將式(5)代入式(1)，則為負壓型液相運動方程。文獻[14]給定下邊界為水頭邊界條件，即z=0時，負壓水頭為h0；上邊界為流量邊界條件，即z=H，比流量為q，建議了如公式(6)的解析解：

(6)

當采用負壓水頭h為變量時,公式(6)中μa-μw的絕對值與負壓水頭絕對值大小相等，將公式(4)及公式(6)化簡合并得：

(7)

公式(7)即為包氣帶一維液相運動方程解析解，該解析解表明土層中有效飽和度剖面受比流量、下邊界負壓水頭、孔隙分布特征和滲透性控制。

若下邊界為潛水位時，則z=0,h0=0,有：

(8)

2.2 均質(zhì)土蒸發(fā)鋒面深度

若規(guī)定比流量q在入滲時為負，蒸發(fā)時為正，當比流量q小于臨界比流量qc時，液態(tài)水飽和度剖面在整個區(qū)域內(nèi)延伸。當比流量q大于臨界比流量qc時，自地表以下至蒸發(fā)鋒面處土孔隙內(nèi)流體為氣體，液態(tài)水僅存在于潛水面至蒸發(fā)鋒面區(qū)間內(nèi)。定義蒸發(fā)鋒面處的有效飽和度為臨界有效飽和度Sec。參考Hayek[3]建議采用在地表處土中水有效飽和度為零時對應(yīng)的蒸發(fā)率為臨界比流量，即:

z=H,Se=0

(9)

將公式(9)代入公式(8)，則有：

(10)

當q>qc時，只需要將臨界有效飽和度Sec代入公式(7)即可推導(dǎo)得蒸發(fā)鋒面深度(蒸發(fā)鋒面距離地下水位的高度)為：

(11)

公式(11)表明蒸發(fā)鋒面深度除受比流量、孔隙分布特征和滲透性影響外，還與臨界有效飽和度Sec有關(guān)。

3 成層土水力連續(xù)區(qū)有效飽和度剖面及蒸發(fā)鋒面深度的解答

均質(zhì)土的分析常用于實驗室研究,成層土常見于實際工程中。根據(jù)Darcy定律和質(zhì)量守恒定律推導(dǎo)包氣帶中成層土液相運動方程的解析解以及相對應(yīng)的蒸發(fā)鋒面深度模型。

3.1 成層土水力連續(xù)區(qū)液相運動方程解析解

現(xiàn)討論兩種滲透性不同的土構(gòu)成非飽和成層土[17]的情形，如圖1所示。令土層的總厚度為L,土層1的厚度為z1,土層2的厚度為L-z1。土層1和土層2的界面位置在z=z1處，每層土為均勻線彈性均質(zhì)土。根據(jù)Darcy定律和質(zhì)量守恒定律，非均質(zhì)界面處土層1和土層2應(yīng)滿足負壓水頭和比流量分別相等。即有：

h1z1=h2z1

(12)

q1z1=q2z1

(13)

式中:h1z1、h2z1分別為非均質(zhì)界面土層1的負壓水頭和土層2的負壓水頭；q1z1、q2z1分別為非均質(zhì)界面處土層1的流量和土層2的流量。

圖1 成層土剖面示意圖

土層1底面為地下水位，故利用式(7),可得土層1中有效飽和度剖面為：

參考譯文:Konka Group has been granted“National Advanced Quality and Benefit Enterprise”and“National Customer Satisfied Enterprise”for several years．

(14)

當z=z1時，根據(jù)負壓水頭和有效飽和度的關(guān)系，可知：

=h2z1

(15)

故對于土層2，其有效飽和度剖面為：

(16)

3.2 成層土蒸發(fā)鋒面深度

對于成層土而言，在穩(wěn)態(tài)蒸發(fā)條件下，蒸發(fā)鋒面深度受比流量影響較大，參考Hayek[3]的建議，本文為便于處理,取土層2頂部有效飽和度為零所對應(yīng)的比流量為土層2的臨界比流量qc2，即z=L,se2=0，代入公式(16)可得土層2的臨界比流量應(yīng)滿足：

qc2=

(17)

定義蒸發(fā)鋒面位于土層非均質(zhì)界面處所對應(yīng)的比流量為臨界比流量qcl。當蒸發(fā)鋒面深度zDF=z1時，即土層2在非均質(zhì)界面處的有效飽和度Se2,z1為土層2的臨界有效飽和度Sec2。將該條件代入公式(16)，可得非均質(zhì)界面處的臨界流量qcl為：

(18)

定義成層土中各土層按所計算土層頂部有效飽和度為臨界飽和度計算臨界比流量，則土層1的臨界比流量應(yīng)滿足當z=z1時Se1=Sec1，代入公式(14)可得：

(19)

當qc2≤q≤qc1時，蒸發(fā)鋒面在土層2內(nèi)，則由公式(16)可得蒸發(fā)鋒面深度為：

(20)

當q≥qc1時，蒸發(fā)鋒面可能在土層1內(nèi)時，若比流量q

(21)

此時，在土層1中0≤z≤zDF內(nèi)有液態(tài)水運動，并服從達西定律。

4 水力參數(shù)反演

解析解重要的一種應(yīng)用是反演水力參數(shù)。水力參數(shù)反演的難易程度取決于反演參數(shù)的個數(shù)。解析式(8)中涉及的水力參數(shù)包括Gardner模型參數(shù)αg和飽和滲透系數(shù)Ks。通過假設(shè)已知土層底面為潛水面來反演上述2個參數(shù)，比流量q、負壓水頭h1、z1處的有效飽和度Se1、z2處的有效飽和度Se2均為已測定值，具體反演過程如下：

將(z1,Se1),(z1,h1)代入公式(4)及公式(8)，化簡可得：

(22)

將(z2,Se2)代入公式(8)求得αg的表達式如下：

(23)

將公式(22)代入公式(23)，整理可得：

(24)

公式(24)是關(guān)于αg的超越方程，需利用迭代求根法求解。

5 結(jié)果與討論

5.1 比流量q和模型參數(shù)αg對有效飽和度剖面的影響

圖2所示為厚度H=10 m土層在不同比流量q下有效飽和度剖面。主要參數(shù)依次為：αg=0.5 m-1,Ks=3×10-7m/s,qc=2.04×10-9m/s，q1=-6×10-8m/s，q2=-6×10-9m/s，q3=1×10-9m/s，q4=6×10-9m/s。根據(jù)圖2可知，當比流量小于臨界比流量時，地表處的有效飽和度隨著比流量的增大而減小，且不同比流量下有效飽和度剖面之間的差異隨距離地下水位高度的減小而減弱。當比流量大于臨界比流量時，土層中存在蒸發(fā)鋒面。土中飽和度剖面以下所圍面積反映了土層中液態(tài)水質(zhì)量分數(shù)相對大小，隨著比流量的增大，土層中液態(tài)水質(zhì)量分數(shù)減小。

圖2 不同比流量下有效飽和度-深度關(guān)系曲線

圖3為厚度H=10 m土層在不同模型參數(shù)αg下有效飽和度剖面。地下水位于土層地面以下10 m處。主要參數(shù)為：αg取值分別為0.05、0.5、1、2 m-1,q=-3×10-8m/s，Ks=3×10-5m/s,由圖3可知αg對有效飽和度剖面影響較大，相同深度處的有效飽和度隨αg值得增大而減小。當αg值較大(如αg=2 m-1)時，地表下部土層一定區(qū)域內(nèi)的有效飽和度近似等于地表處的有效飽和度，在該區(qū)域以下有效飽和度會驟增，直至地下水位處變?yōu)?。此時液相運動主要受重力控制。當αg值較小(如αg=0.05 m-1)時，毛細作用發(fā)揮支配作用，有效飽和度趨于線性變化。

圖3 不同αg下有效飽和度-深度關(guān)系曲線

5.2 成層土中有效飽和度剖面的變化特征

為分析穩(wěn)態(tài)流下成層土有效飽和度隨深度的變化，設(shè)定土層A和土層B的厚度均為5 m。土層A的主要參數(shù)為：αg=1 m-1，Ks=3×10-5m/s;土層B的參數(shù)為:αg=0.5 m-1，Ks=3×10-7m/s。土層A具有比土層B較高的孔隙連通狀態(tài)和透水性。

5.2.1 入滲條件下成層土有效飽和度剖面 圖4為穩(wěn)態(tài)入滲條件下成層土內(nèi)有效飽和度剖面。其中工況1為土層B在上部，土層A在下部；工況2為土層A在上部，土層B在下部。比流量q=-3×10-8m/s。

工況1情況下，上部土層B的有效飽和度小于相同深度處均質(zhì)土B的有效飽和度；工況2情況下，上部土層A的有效飽和度大于相同深度處均質(zhì)土A的有效飽和度。而對于下部土層，由于下邊界水頭與均質(zhì)土下邊界條件相同，無論是工況1還是工況2均與對應(yīng)均質(zhì)土層的有效飽和度-深度關(guān)系曲線重合。

圖4表明穩(wěn)態(tài)入滲下，當下部土層透水性較好時，上部土層中的液態(tài)水向下層流動，液相運動增強，其有效飽和度與均質(zhì)土相比要小。反之，當下部土層透水性較差時，上部土層中的液態(tài)水向下層流動受到阻滯，液相運動減弱，其有效飽和度與均質(zhì)土相比要大。

圖4 入滲條件下成層土有效飽和度-深度關(guān)系曲線

5.2.2 蒸發(fā)條件下成層土有效飽和度剖面 圖5為穩(wěn)態(tài)蒸發(fā)條件下成層土的有效飽和度剖面。其中，工況3為土層B在上部，土層A在下部；工況4為土層A在上部，土層B在下部。比流量q=1×10-8m/s。由圖5可知，在蒸發(fā)條件下工況3的上部土層B內(nèi)有效飽和度大于相同深度處均質(zhì)土B的有效飽和度；工況4的上部土層A的有效飽和度小于相同深度處均值土A的有效飽和度。下部土層中，由于下邊界水頭與均質(zhì)土下邊界條件相同，工況3、工況4均與對應(yīng)的均質(zhì)土層的有效飽和度-深度關(guān)系曲線重合。

對比上部土層的有效飽和度，可以發(fā)現(xiàn)在穩(wěn)態(tài)蒸發(fā)下，若下部土層透水性較好，上部土層中液相運動增強會儲存較多的液態(tài)水，同時蒸發(fā)鋒面距地表更淺；反之，當下部土層透水性較差時，上部土層中液相運動受阻減弱，儲存液態(tài)水較少，蒸發(fā)鋒面距地表更深。

圖5 蒸發(fā)條件下成層土有效飽和度-深度關(guān)系曲線

5.3 參數(shù)反演分析

5.3.1 參數(shù)反演模型的可靠性 假設(shè)土層底面為地下水位，給定q=-3×10-8m/s,αg=0.5 m-1，Ks=3×10-7m/s。選擇兩個觀測點Z1=8 m,Z2=4 m,利用公式(8)獲得Z1處有效飽和度精確解為0.116 5，將Z1處精確解代入公式(5)中，可得負壓水頭的精確解為-4.299 73；Z2處的有效飽和度精確解為0.221 8。將上述數(shù)據(jù)代入公式(24)中求得αg值為0.5 m-1,再代入公式(22)求得Ks值為3×10-7m/s。結(jié)果表明參數(shù)反演模型的正確。

6 結(jié) 論

本文利用相對滲透系數(shù)模型中的Gardner模型分別推導(dǎo)了均質(zhì)非飽和土和成層非飽和土中一維液相運動方程解析解，獲得了蒸發(fā)鋒面深度的計算表達式，所得結(jié)論具體如下：

(1)包氣帶中均質(zhì)土有效飽和度是深度的函數(shù)，受比流量、下邊界負壓水頭、Gardner模型參數(shù)等參數(shù)控制。相同深度處的有效飽和度隨著比流量和Gardner模型參數(shù)的增大而減小。

(2)對于給定比流量和土層參數(shù)的非飽和成層土中，土層內(nèi)的有效飽和度剖面主要受下邊界負壓水頭控制。在穩(wěn)態(tài)入滲或穩(wěn)態(tài)蒸發(fā)條件下，下部土層的有效飽和度剖面與相應(yīng)均質(zhì)土的有效飽和度剖面重合；上部土層的有效飽和度受下部土層透水性的影響。

(3)由于蒸發(fā)鋒面的存在，包氣帶中均質(zhì)土和成層土液相運動方程的解析解的適用范圍在臨界有效飽和度與飽和有效飽和度之間。

(4)測試了參數(shù)反演模型的可靠性，分析測量數(shù)據(jù)誤差對參數(shù)反演結(jié)果影響。結(jié)果表明所建議的參數(shù)反演表達式是可靠的，參數(shù)反演的結(jié)果受測量誤差的影響較小。