兩種快速測(cè)定非飽和土 水力學(xué)參數(shù)方法的對(duì)比分析

伊盼盼 ，韋昌富，魏厚振，曹華峰，徐文強(qiáng)

（1.中國科學(xué)院武漢巖土力學(xué)研究所 巖土力學(xué)與工程國家重點(diǎn)實(shí)驗(yàn)室，武漢 430071；2.武昌理工學(xué)院，武漢 430223）

1 引 言

非飽和土水力學(xué)參數(shù)包括土-水特征曲線和滲透函數(shù)。土-水特征曲線（SWCC）表示非飽和土的基質(zhì)吸力與含水率（如飽和度或體積含水率）之間的關(guān)系。滲透函數(shù)是滲透系數(shù)與含水率之間的關(guān)系。這兩大水力學(xué)參數(shù)在分析非飽和土變形與強(qiáng)度、降雨入滲與邊坡穩(wěn)定、污染物地下遷移等問題中扮演著重要角色[1－2]。

傳統(tǒng)的非飽和土水力學(xué)參數(shù)的測(cè)定方法均是在試樣處于平衡態(tài)下測(cè)得的，耗時(shí)較長(zhǎng)。另外，由于加載時(shí)間較長(zhǎng)，高壓氣體會(huì)通過水的流動(dòng)在陶土板背面析出，影響測(cè)試精度。因此，如何能夠快速、有效地確定非飽和土水力學(xué)參數(shù)成為需要解決的關(guān)鍵問題。

通過流動(dòng)試驗(yàn)并利用數(shù)值反算方法對(duì)結(jié)果進(jìn)行分析能夠快速得到非飽和土的水力學(xué)參數(shù)。Gardner[3]首次提出流動(dòng)試驗(yàn)法，該方法通過求解非線性流動(dòng)方程得到非飽和土水力學(xué)參數(shù)。Kool 等[4]提出一步流動(dòng)試驗(yàn)方法，結(jié)合數(shù)值反算方法得到非飽和土水力學(xué)參數(shù)。盡管一步流動(dòng)試驗(yàn)?zāi)軐?duì)試驗(yàn)結(jié)果進(jìn)行很好地?cái)M合，但容易出現(xiàn)參數(shù)不惟一的現(xiàn)象。van Dam[5－6]等指出若能同時(shí)測(cè)得試樣中的孔隙水壓力，則可以降低參數(shù)不惟一的問題。另外，Eching等[7－9]同時(shí)開展了一步流動(dòng)和多步流動(dòng)試驗(yàn)，結(jié)果顯示多步流動(dòng)試驗(yàn)?zāi)芨玫慕鉀Q參數(shù)不惟一的問題。

文中提出能夠快速測(cè)定非飽和土水力學(xué)參數(shù)的一步流動(dòng)和動(dòng)態(tài)多步流動(dòng)方法，并對(duì)這兩種方法進(jìn)行了對(duì)比分析。其中，一步流動(dòng)方法是利用HYDRUS-1D 一維水分運(yùn)移模型擬合試樣一步脫濕流動(dòng)試驗(yàn)下溢出水量隨時(shí)間變化的關(guān)系曲線，進(jìn)而反算得出非飽和土的土-水特征曲線和滲透函數(shù)。該方法能夠降低參數(shù)的不惟一性，同時(shí)考慮了陶土板材料參數(shù)對(duì)試驗(yàn)溢出水隨時(shí)間變化規(guī)律的影響。動(dòng)態(tài)多步流動(dòng)方法與上述多步流動(dòng)試驗(yàn)不同，它是基于Wei 和Muraleetharan[10]、Wei 和Dewoolkar[11]提出的能夠描述靜態(tài)土-水特征曲線和動(dòng)態(tài)土-水特征曲線間關(guān)系的動(dòng)態(tài)模型，開展多步流動(dòng)試驗(yàn)測(cè)定非飽和土的水力學(xué)參數(shù)。在此動(dòng)態(tài)多步流動(dòng)試驗(yàn)過程中不需要土樣處于水力學(xué)平衡狀態(tài)，因而節(jié)省了試驗(yàn)所需時(shí)間。該試驗(yàn)中試樣在每一吸力步下都達(dá)到平衡。

2 兩種方法的原理

2.1 動(dòng)態(tài)多步流動(dòng)方法

2.1.1 動(dòng)態(tài)（非平衡態(tài)）土-水特征關(guān)系

考慮一圓柱形非飽和土樣，其兩端分別與氣源和水源相連接。在初始時(shí)，土樣處于靜力平衡狀態(tài)，其飽和度為，吸力為=，這里( Sr)是平衡態(tài)的土-水特征函數(shù)。此時(shí)，若改變?cè)嚇觾啥说膲毫Σ?，試樣中的水?huì)從一端流向另外一端，故試樣的飽和度也會(huì)隨之發(fā)生變化。Wei 和 Muralee- tharan[10]與Wei 和Dewoolkar[11]提出基質(zhì)吸力的變化與飽和度之間存在以下關(guān)系：

式中：基質(zhì)吸力 pc= pa- pw，這里 pa為孔隙氣壓力， pw為孔隙水壓力；μ 為材料參數(shù)；θ 為體積含水率。

如果不考慮試樣的體積變化，式（1）可以轉(zhuǎn)化為

式中： nμ μ′= ，n 為孔隙率。

若吸力增加一小量 cpδ ，則飽和度的變化 rSδ與 cpδ 的關(guān)系為

式（3）類似于線性、黏彈性固體本構(gòu)關(guān)系Kelvin模型，可以用圖1 來表示。由此可見，當(dāng)吸力改變時(shí)，飽和度不是突然變化的，而是隨時(shí)間慢慢地改變，這一過程與實(shí)際是相吻合的。

圖1 線性黏彈性Kelvin模型 Fig.1 Linear viscoelastic model of Kelvin

假設(shè)擾動(dòng)（吸力）是突然增加的，故有

式中： ( )H t 為Heaviside 函數(shù)。

那么求解式（3），得

如果吸力是多步施加的，即

同理，求解式（3）可以得到多步流動(dòng)下試樣飽和度隨時(shí)間的演化方程為

式中：kC 和分別為吸力為時(shí)的穩(wěn)態(tài)土-水特征曲線斜率（容水率）和擾動(dòng)后孔隙水流動(dòng)的特征時(shí)間。

2.1.2 平衡條件下土-水特征曲線和滲透系數(shù)的確定

由式（7）可以看出，式中主要由C 和cτ 兩大參數(shù)。如果這兩個(gè)參數(shù)值已知，那么試樣的飽和度隨時(shí)間的演化曲線也就明確了，因此，確定這兩大參數(shù)值非常重要。這里列出兩種求解的方法。

方法1 為兩點(diǎn)法，考慮某時(shí)間段 ti-1＜ t ＜ ti中的兩個(gè)時(shí)刻和，由式（7）可以得到一個(gè)非線性方程組，采用迭代法可以求出每級(jí)吸力步對(duì)應(yīng)的C 和τc值。

方法2 為最小二乘法。根據(jù)實(shí)測(cè)的飽和度隨時(shí)間的演化曲線，在上面取足夠多具有代表性的點(diǎn)，采用最小二乘法對(duì)實(shí)測(cè)數(shù)據(jù)進(jìn)行擬合，求出每級(jí)吸力步對(duì)應(yīng)的C 和cτ 值。

確定了各級(jí)吸力作用下的容水率C 和特征時(shí)間τc，就可以得到在基質(zhì)吸力作用下試樣處于平衡態(tài)時(shí)對(duì)應(yīng)的飽和度：

對(duì)于一定尺寸的土樣，孔隙水壓力耗散時(shí)間（特征時(shí)間）cτ 與土樣的滲透性有關(guān)，通過對(duì)孔隙水的滲流過程分析，可以建立特征時(shí)間cτ 與滲透系數(shù)的關(guān)系。其過程如圖2 所示。由圖可知，初始時(shí)土樣處于平衡狀態(tài)，在某個(gè)時(shí)刻，土樣底部的孔隙水壓力降低了 wpδ ，而其他條件不變。為了描述這一過程，假定在整個(gè)過程中孔隙氣體壓力ap 保持不變，同時(shí)假定土顆粒和水是不可壓縮的（即質(zhì)量密度sρ 和wρ 為常數(shù)），且不考慮變形的影響（孔隙率為常數(shù)）。

圖2 土中Darcy流動(dòng)示意圖 Fig.2 Darcy flow diagram of the soil

在這些前提下，孔隙水質(zhì)量平衡方程為

式中：w

n 為水相孔隙率。

式（9）可變形為

孔隙水的流動(dòng)由以下Darcy 類方程控制：

利用土-水特征關(guān)系rc( )S p 及C 的定義，將式（11）代入式（10）可得：

亦可以寫成：

這就是在上述試驗(yàn)條件下土樣孔隙水壓力耗散過程的控制方程，式中vc 為孔隙水壓力耗散系數(shù)。

設(shè)試樣的長(zhǎng)度為L(zhǎng)，那么由式（13）得出Darcy流動(dòng)的特征時(shí)間cτ 為

2.2 一步流動(dòng)方法

一步流動(dòng)方法的原理為：通過對(duì)飽和試樣施加一步吸力，測(cè)得試樣溢出水量隨時(shí)間的變化曲線。利用一維水分運(yùn)移過程分析軟件HYDRUS-1D 對(duì)溢出水隨時(shí)間流量曲線進(jìn)行擬合，同時(shí)分析反算出一些難以實(shí)測(cè)的水力學(xué)參數(shù)，進(jìn)而得到非飽和土-水特征曲線和滲透函數(shù)曲線。

土樣中水分變化可由Richards 方程描述，即

式中：h 為壓力水頭； ( )K h 為滲透函數(shù)；t 為時(shí)間；z 為土樣深度。

試樣的初邊值條件表示為

對(duì)試樣施加一步吸力P 之前，先施加一小吸力0P （略大于試樣進(jìn)氣值），使試樣由飽和狀態(tài)變?yōu)榉秋柡蜖顟B(tài)。這里00/h P gρ= ， 0z= 是試樣上部邊界，z L= 為試樣下部邊界， ( )h t 為試樣下部邊界的水頭。

土-水特征曲線模型選用比較常用的van Genu- chten[12]模型： ）

將van Genuchten 模型與Mualem[13]模型相結(jié)合可得到滲透函數(shù)的表達(dá)式為

式中：sθ 為飽和含水率；rθ 為殘余含水率；α 、β 、γ 為經(jīng)驗(yàn)參數(shù)；sK 為飽和滲透系數(shù)；eS 為有效含水率；l 為孔隙結(jié)構(gòu)參數(shù)。

3 試驗(yàn)設(shè)備和試驗(yàn)方法

3.1 試驗(yàn)設(shè)備

為了開展前面提出的多步流動(dòng)試驗(yàn)，采用聯(lián)合測(cè)試系統(tǒng)，其示意圖如圖3 所示。

圖3 聯(lián)合測(cè)試系統(tǒng)示意圖 Fig.3 Sketch of combined testing system

該系統(tǒng)主要由壓力控制系統(tǒng)、壓力室、儲(chǔ)水容器、稱量系統(tǒng)、數(shù)據(jù)采集系統(tǒng)等主要部分構(gòu)成。

壓力控制系統(tǒng)由高壓氮?dú)馄?、減壓閥、壓力調(diào)節(jié)泵等部分構(gòu)成。減壓閥把氮?dú)馄坷锏母邏簹怏w變成低壓氣體，起粗調(diào)作用。壓力調(diào)節(jié)泵能精確控制給壓力室輸入的氣壓力值，起到細(xì)調(diào)作用。

壓力室主要由裝樣容器和上、下蓋板等3 部分組成，其中下蓋板裝有高進(jìn)氣值陶土板。

儲(chǔ)水容器的用途有兩個(gè)，一是頂部和真空泵相連，可以用來對(duì)試樣進(jìn)行飽和；二是用于對(duì)系統(tǒng)管路的飽和。

稱量系統(tǒng)主要由天平及上方開有一小孔的盛水杯組成，開一小孔有利于防止水分揮發(fā)。

數(shù)據(jù)采集系統(tǒng)是通過LabVIEW 軟件來完成數(shù)據(jù)采集功能，它有兩個(gè)數(shù)據(jù)通道分別與天平和壓力調(diào)節(jié)泵相連，可以采集試樣的排出水量和對(duì)試樣施加的吸力隨時(shí)間變化的數(shù)據(jù)。采集的時(shí)間間隔可以從幾秒到幾小時(shí)不等，可根據(jù)試驗(yàn)需要任意設(shè)置。

3.2 試驗(yàn)方法

試驗(yàn)所用土樣取自黃河三角洲，其基本物理特性指標(biāo)如表1 所示。

表1 土樣基本物理性質(zhì) Table 1 Physical properties of the soil sample

取該土樣以初始含水率為9%，孔隙率以42%為控制目標(biāo)，制備一個(gè)直徑為5.4 cm，高為3 cm 的圓柱形試樣。隨后，將制備好試樣采用抽氣飽和法進(jìn)行飽和。先對(duì)飽和試樣開展多步流動(dòng)試驗(yàn)，具體過程如下：

在 t0時(shí)刻，對(duì)試樣施加第1 級(jí)吸力增量，然后在吸力保持不變的條件下，直至 t1時(shí)刻（此時(shí)試樣并沒有達(dá)到平衡狀態(tài)）。接著，在 t1時(shí)刻，施加第2 級(jí)吸力增量，并保持吸力不變至 t2時(shí)刻。以此類推，施加第N 級(jí)吸力增量，直至 tn時(shí)刻。根據(jù)這組多步流動(dòng)試驗(yàn)數(shù)據(jù)，利用前面的理論模型及求解辦法，可以推算出試樣脫濕條件下平衡態(tài)的土-水特征曲線。

由于文中采用的試樣在試驗(yàn)過程中其體變可忽略不計(jì)，且經(jīng)歷了脫濕、吸濕之后試樣的結(jié)構(gòu)基本趨于穩(wěn)定，故在多步流動(dòng)試驗(yàn)后把試樣重新飽和，然后對(duì)其開展一步流動(dòng)試驗(yàn)，初始飽和度的影響可以忽略，其過程如下：對(duì)該飽和試樣施加一小吸力步（略大于試樣的進(jìn)氣值），使之由飽和狀態(tài)變?yōu)榉秋柡蜖顟B(tài)。當(dāng)流動(dòng)停滯（穩(wěn)定狀態(tài)）時(shí)，對(duì)試樣施加一較大吸力步，直至溢出水量穩(wěn)定。這樣就可以得到脫濕段溢出水量隨時(shí)間變化的關(guān)系曲線。利用一維水分運(yùn)移過程分析軟件HYDRUS-1D 分別對(duì)脫濕段和吸濕段溢出水隨時(shí)間的流量曲線進(jìn)行擬合，反算得出一些難以實(shí)測(cè)的水力學(xué)參數(shù)，進(jìn)而得到非飽和土干、濕土-水特征曲線和滲透函數(shù)曲線。

4 試驗(yàn)結(jié)果分析

4.1 動(dòng)態(tài)多步流動(dòng)試驗(yàn)結(jié)果

由于每級(jí)吸力下加載的時(shí)間不同，測(cè)得的試樣飽和度隨時(shí)間的演化關(guān)系也不同，本文分別進(jìn)行了兩組不同時(shí)間步的流動(dòng)試驗(yàn)。每組流動(dòng)試驗(yàn)施加的吸力步及各吸力步的加載時(shí)間如表2 所示。為了進(jìn)行對(duì)比分析，試樣在各吸力步下達(dá)到平衡的時(shí)間也在表2 中列出。

由表可以看出，這兩組流動(dòng)試驗(yàn)在各吸力下的加載時(shí)間遠(yuǎn)小于試樣達(dá)到平衡所需的時(shí)間，第2 組流動(dòng)試驗(yàn)的加載時(shí)間比第1 組稍長(zhǎng)。另外，這兩組流動(dòng)試驗(yàn)在初級(jí)吸力下都達(dá)到了平衡，原因是試樣由飽和區(qū)向非飽和區(qū)過渡時(shí)，土-水特征曲線的斜率即C 值變化較大，故先對(duì)試樣施加一小吸力（略大于進(jìn)氣值），待達(dá)到平衡時(shí)，試樣處于非飽和狀態(tài)，以該狀態(tài)為初始狀態(tài)開展多步流動(dòng)試驗(yàn)?zāi)苡行У臏p少試驗(yàn)誤差。

這兩組流動(dòng)試驗(yàn)得到的試樣飽和度隨時(shí)間的演化曲線如圖4 中數(shù)據(jù)點(diǎn)所示，利用前面提出的動(dòng)態(tài)模型，并采用求解方法2，分別對(duì)這兩組試驗(yàn)曲線進(jìn)行最小二乘擬合，其擬合結(jié)果如圖4 中實(shí)線所示，最小二乘擬合的精度如圖5 所示。

從以上圖4、5 可以看出，這兩組流動(dòng)試驗(yàn)的最小二乘擬合精度都比較高，進(jìn)而說明了由最小二乘擬合得出的各級(jí)吸力步下C、τ 值的可靠性。

4.2 一步流動(dòng)試驗(yàn)結(jié)果

先對(duì)飽和試樣施加一較小吸力1 kPa，有少量水溢出，這時(shí)試樣從飽和狀態(tài)進(jìn)入非飽和狀態(tài)，然后增加吸力至80 kPa，進(jìn)行脫濕試驗(yàn)。待溢出水量恒定后，根據(jù)試樣初始飽和質(zhì)量和溢出水的質(zhì)量可以得到體積含水率隨時(shí)間變化的關(guān)系。

根據(jù)該關(guān)系曲線、初邊值條件、陶土板材料參數(shù)運(yùn)用HYDRUS-1D 一維水分運(yùn)移模型進(jìn)行反算，擬合出試樣脫濕段溢出水量隨時(shí)間的演化曲線。實(shí)測(cè)數(shù)據(jù)曲線和擬合曲線如圖6 所示。

圖4 兩組流動(dòng)試驗(yàn)中實(shí)測(cè) S r-t 曲線與擬合結(jié)果的對(duì)比 Fig.4 Comparison between measured and fitted curves of S r-t during the two outflow tests

圖5 兩組流動(dòng)試驗(yàn)中最小二乘的擬合精度 Fig.5 Accuracy of least square fitting during the two outflow tests

圖6 脫濕段θ -t 曲線模型預(yù)測(cè)值與實(shí)測(cè)值比較 Fig.6 Comparison between estimated and measured values of θ -t curves during drying

從圖中可以看出，HYDRUS-1D 內(nèi)部模型能夠很好地?cái)M合實(shí)測(cè)數(shù)據(jù)，其擬合參數(shù)如表3 所示。

表3 脫濕條件下的擬合參數(shù) Table 3 Fitting parameters for drying

4.3 試驗(yàn)結(jié)果對(duì)比分析

對(duì)于多步流動(dòng)試驗(yàn)，由于各級(jí)吸力步下的C、τ 值已知，根據(jù)前面的理論模型，可得到平衡狀態(tài)下試樣的土-水特征曲線和滲透函數(shù)曲線；對(duì)于一步流動(dòng)試驗(yàn)，根據(jù)擬合參數(shù)可反算得出平衡狀態(tài)下試樣的土-水特征曲線和滲透函數(shù)曲線。由這兩種試驗(yàn)方法得到的試樣在平衡態(tài)下的土-水特征曲線和實(shí)測(cè)值的對(duì)比結(jié)果如圖7 所示。

從圖可以看出，一步流動(dòng)試驗(yàn)得到的試樣平衡態(tài)土-水特征曲線與實(shí)測(cè)值有些偏差，而多步流動(dòng)試驗(yàn)結(jié)果則與實(shí)測(cè)值接近，第2 組流動(dòng)試驗(yàn)的加載時(shí)間較長(zhǎng)，得到的土-水特征曲線與實(shí)測(cè)值基本吻合。

由這兩種試驗(yàn)方法得到的試樣滲透函數(shù)曲線，如圖8 所示。

圖7 兩種流動(dòng)試驗(yàn)及實(shí)測(cè)得到的SWCC對(duì)比 Fig.7 Comparison of SWCC between the measured and the two outflow tests

圖8 兩種流動(dòng)試驗(yàn)得到的滲透函數(shù)曲線對(duì)比 Fig.8 Comparison curves of hydraulic permeability function between the two outflow tests

由圖可以看出，這兩種試驗(yàn)方法得到的滲透函數(shù)接近。對(duì)比多步流動(dòng)結(jié)果，一步流動(dòng)試驗(yàn)反算得到的滲透系數(shù)在試樣體積含水率高時(shí)偏高，在體積含水率低時(shí)偏低。這是由于一步流動(dòng)方法是根據(jù)溢出水量隨時(shí)間的變化關(guān)系反算得到，受邊界條件和陶土板材料參數(shù)的制約。

5 結(jié) 論

傳統(tǒng)的測(cè)試非飽和土力學(xué)參數(shù)的方法，需要很長(zhǎng)時(shí)間，并且難以克服由于加載時(shí)間較長(zhǎng)高壓氣體在陶土板背面析出的問題。因此，快速測(cè)定非飽和水力學(xué)參數(shù)的方法研究顯得尤為重要，通過研究可得如下結(jié)論：

（1）一步流動(dòng)方法和動(dòng)態(tài)多步流動(dòng)方法均能在短時(shí)間內(nèi)測(cè)定非飽和土的水力學(xué)參數(shù)，克服了傳統(tǒng)測(cè)試方法耗時(shí)較長(zhǎng)的問題，并且有效地避免了高壓氣體通過水的流動(dòng)在陶土板背面析出的問題，提高了測(cè)定非飽和土水力學(xué)參數(shù)的速度和準(zhǔn)確性。

（2）動(dòng)態(tài)多步流動(dòng)方法得到試樣平衡態(tài)的土-水特征曲線與實(shí)測(cè)值基本吻合，無需等到試樣達(dá)到平衡狀態(tài)后再進(jìn)行測(cè)定，大大縮短了測(cè)試時(shí)間；一步流動(dòng)方法僅對(duì)試樣施加一步吸力，比動(dòng)態(tài)多步流動(dòng)方法所用的時(shí)間短，但得出的試樣平衡態(tài)土-水特征曲線和實(shí)測(cè)數(shù)值有些偏差。因此，通過對(duì)比分析，動(dòng)態(tài)多步流動(dòng)方法的結(jié)果要優(yōu)于一步流動(dòng)試驗(yàn)方法。

[1] FREDLUND D G, XING A Q, FREDLUND M D, et al. The relationship of the unsaturated soil shear strength function to the soil-water characteristic curve[J]. Cana- dian Geotechnical Journal, 1996, 32(3): 440－448.

[2] VANAPALLI S K, FREDLUND D G, PUFAHL D E, et al. Model for the prediction of shear strength with respect to soil suction[J]. Canadian Geotechnical Journal, 1996, 33(3): 379－392.

[3] GARDNER W R. Calculation of capillary conductivity from pressure plate outflow data[J]. Soil Science Society of America Journal, 1956, 20: 317－320.

[4] KOOL J B, PARKER J C, VAN GENUCHTEN M TH. Determining soil hydraulic properties from one-step outflow experiments by parameter estimation: I. Theory and numerical studies[J]. Soil Science Society of America Journal, 1985, 49: 1348－1354.

[5] VAN DAM J C, STRICKER J N M, DROOGERS P. Evaluation of the inverse method for determining soil hydraulic functions from one-step outflow experiments[J]. Soil Science Society of America Journal, 1992, 56: 1042－1050.

[6] TOORMAN A F, WIERENGA P J, HILLS R G. Parameter estimation of hydraulic properties from one- step outflow data[J]. Water Resources Research, 1992, 28(11): 3021－3028.

[7] ECHING S O, HOPMANS J W. Optimization of hydraulic functions from transient outflow and soil-water pressure data[J]. Soil Science Society of America Journal, 1993, 57(5): 1167－1175.

[8] ECHING S O, HOPMANS J W, WENDROTH O. Unsaturated hydraulic conductivity from transient multistep outflow and soil-water pressure data[J]. Soil Science Society of America Journal, 1994, 58: 687－695.

[9] VAN DAM J C, STRICKER J N M, DROOGERS P. Inverse method to determine soil hydraulic functions from multistep outflow experiments[J]. Soil Science Society of America Journal, 1994, 58: 647－652.

[10] WEI C F, DEWOLKAR M M. Formulation of capillary hysteresis with internal state variables[J]. Water Resources Research, 2006, 42(7): 1－16.

[11] WEI C F, MURALEETHARAN K K. Linear viscoelastic behavior of porous media with non-uniform saturation[J]. International Journal of Engineering Science, 2007, 45(2－8): 698－715.

[12] VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44: 892－898.

[13] MUALEM Y. A new model for predicting the hydraulic conductivity of unsaturated porous media[J]. Water Resources Research, 1976, 12(3): 513－522.