水-氣二相流本構(gòu)模型參數(shù)的反演識(shí)別

孫冬梅 ，馮 平，張明進(jìn)，Semprich S

(1. 天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；2. 交通運(yùn)輸部天津水運(yùn)工程科學(xué)研究所水工構(gòu)造物檢測(cè)、診斷與加固技術(shù)實(shí)驗(yàn)室，天津 300456；3. 奧地利格拉茨技術(shù)大學(xué)土力學(xué)及基礎(chǔ)工程研究所，格拉茨 A-8010)

目前，多相流問(wèn)題的研究在水文地質(zhì)、石油工程及環(huán)境工程等領(lǐng)域迅速發(fā)展起來(lái)，多相流模型的進(jìn)一步發(fā)展對(duì)地下水資源的模擬和評(píng)價(jià)、對(duì)已污染地下水土的治理與修復(fù)，以及地下水環(huán)境災(zāi)害的防治等問(wèn)題具有重要意義[1-2]．在多相流數(shù)值模擬中，需要首先確定與各相流體(包括水相、氣相和非水相流體（NAPL)）有關(guān)的本構(gòu)關(guān)系，包括毛細(xì)壓力-飽和度、相對(duì)滲透率-飽和度關(guān)系，表征這些本構(gòu)關(guān)系的模型包括 Brooks-Corey(BC)模型和 Brooks-Corey-Mualem(BCM)模型、van Genuchten(VG)模型和 van Genuchten-Mualem(VGM)模型等[3-5]，這些模型中包含一組待定的與多孔介質(zhì)孔隙結(jié)構(gòu)特性有關(guān)的參數(shù)．在多相流情況下，同時(shí)確定水相、氣相和非水相流體間的本構(gòu)關(guān)系模型參數(shù)是十分困難的，基于Leverett和Lewis的定性描述，水-氣二相流體的本構(gòu)關(guān)系僅與各自的飽和度有關(guān)，因而目前較為有效的研究方法是首先確定水-氣二相流體的本構(gòu)關(guān)系模型參數(shù)，以此為基礎(chǔ)，推求其他流體對(duì)(NAPL-氣、水-NAPL)間的本構(gòu)關(guān)系模型參數(shù)[6-7]，因此，確定水-氣二相流體的本構(gòu)關(guān)系模型參數(shù)是多相流模型研究中最基本的前提．

確定水-氣二相流體的本構(gòu)關(guān)系模型參數(shù)較為常用的方法是直接測(cè)量法，該方法主要通過(guò)平衡、穩(wěn)態(tài)的試驗(yàn)來(lái)確定本構(gòu)模型參數(shù)，最大的困難是進(jìn)行穩(wěn)態(tài)的試驗(yàn)通常受到初始條件和邊界條件的限制，而且比較費(fèi)時(shí)．近年來(lái)，隨著數(shù)值計(jì)算水平及試驗(yàn)技術(shù)的發(fā)展，建立水-氣二相流數(shù)學(xué)模型及進(jìn)行有效的非穩(wěn)定水-氣二相流試驗(yàn)成為可能，采用反演識(shí)別的方法的優(yōu)勢(shì)越來(lái)越突顯[8-10]，其主要優(yōu)點(diǎn)是數(shù)值反演的水-氣二相流試驗(yàn)可以是瞬態(tài)的、非穩(wěn)定的流動(dòng)過(guò)程，相比平衡、穩(wěn)態(tài)的二相流試驗(yàn)要靈活、省時(shí)，而且試驗(yàn)的邊界條件和初始條件可以自由設(shè)置，另外可以同時(shí)確定一組表征同一土樣的水-氣二相流體的本構(gòu)關(guān)系模型參數(shù)．

采用反演識(shí)別的方法來(lái)確定水-氣二相流本構(gòu)模型參數(shù)的過(guò)程中，首先要進(jìn)行一個(gè)非穩(wěn)定水-氣二相流試驗(yàn)來(lái)獲取試驗(yàn)數(shù)據(jù)，試驗(yàn)的設(shè)計(jì)非常關(guān)鍵，應(yīng)使得反演模型的解具有唯一性和穩(wěn)定性．如 Gardner[11]所設(shè)計(jì)的單步非穩(wěn)定空氣出流試驗(yàn)就未能使得反演模型的解具有唯一性，主要原因是試驗(yàn)數(shù)據(jù)的信息不足；在此基礎(chǔ)上，van Dam 等[12]進(jìn)行了多步的非穩(wěn)定空氣出流試驗(yàn)，并且測(cè)量了試驗(yàn)過(guò)程中的累積出流量數(shù)據(jù)，從而獲得了足夠的試驗(yàn)信息使得反演模型的解具有唯一性；后來(lái)經(jīng)過(guò) Eching和 Hopmans[13]的研究表明，在非穩(wěn)定水-氣二相流試驗(yàn)中測(cè)量孔隙壓力和累積出流量數(shù)據(jù)，能夠使得反演模型的解具有唯一性和穩(wěn)定性．因此確定水-氣二相流本構(gòu)模型參數(shù)的反演識(shí)別方法需要的準(zhǔn)備工作主要包括3部分：設(shè)計(jì)合理的非穩(wěn)定水-氣二相流試驗(yàn)、能夠有效模擬水-氣二相流動(dòng)的水-氣二相流數(shù)值模擬模型以及基于最優(yōu)化算法所建立的本構(gòu)模型參數(shù)的反演識(shí)別模型．可見(jiàn)，進(jìn)行水-氣二相流本構(gòu)模型參數(shù)的反演識(shí)別研究過(guò)程中存在許多困難，該研究是多相流數(shù)值模擬研究進(jìn)一步發(fā)展的瓶頸問(wèn)題．

筆者進(jìn)行了多步的非穩(wěn)定空氣出流試驗(yàn)，測(cè)量了試驗(yàn)過(guò)程中兩點(diǎn)的孔隙氣壓力、孔隙水壓力、氣相飽和度及累積出流量的變化過(guò)程；建立了考慮水相、氣相流動(dòng)及相互溶混的水-氣二相流數(shù)值模擬模型，采用 Marquart-Levenberg最優(yōu)化算法建立了水-氣二相流本構(gòu)關(guān)系模型參數(shù)反演識(shí)別模型，結(jié)合非穩(wěn)定水-氣二相流試驗(yàn)，對(duì)水-氣二相流的本構(gòu)關(guān)系模型(VG和VGM)參數(shù)進(jìn)行反演識(shí)別．

1 非穩(wěn)定水-氣二相流試驗(yàn)

本研究所采用的試驗(yàn)裝置如圖1所示，土樣放置在圓筒狀樹(shù)脂玻璃容器的中部，土樣高度為 0.45,m，其中距離底部 1/3處和頂部 1/3處放置時(shí)間域反射計(jì)TDR(T1～T4)和壓力表(P2～P3)，分別測(cè)量土柱在試驗(yàn)過(guò)程中的體積含水量和土柱中孔隙氣壓力．TDR和壓力表都與數(shù)據(jù)采集器相連，試驗(yàn)過(guò)程數(shù)據(jù)通過(guò)計(jì)算機(jī)輸出．頂部溢流出水量由電子秤測(cè)量其質(zhì)量．底部與供氣裝置相連接，壓力表 P1測(cè)量底部入流的氣體壓力．

圖1 非穩(wěn)定水-氣二相流試驗(yàn)裝置Fig.1 Layout of unsteady water-air two-phase flow test setup

試驗(yàn)方法：首先將土樣在攪拌機(jī)中攪拌均勻，然后將土樣一層一層地放到樹(shù)脂玻璃圓筒容器中，分層壓實(shí)，測(cè)定此時(shí)土樣的干密度和孔隙率．在分層壓實(shí)的過(guò)程中，在相應(yīng)的位置放置 TDR和壓力表．通過(guò)對(duì)土樣進(jìn)行一系列常水頭滲透試驗(yàn)，測(cè)得土樣的固有滲透系數(shù)和飽和水飽和度．利用供水裝置，分級(jí)增大土樣頂部的水頭，使土樣從上而下漸進(jìn)地達(dá)到飽和.最后當(dāng)TDR顯示的含水量達(dá)到穩(wěn)定狀態(tài)時(shí)即達(dá)到飽和，最后一直供水，直到土樣上部的水頭穩(wěn)定不變，此時(shí)，土樣達(dá)到穩(wěn)定飽和狀態(tài)，氣驅(qū)替水的空氣出流試驗(yàn)開(kāi)始進(jìn)行．關(guān)掉頂部的供水管，打開(kāi)底部的供氣管，可調(diào)節(jié)的供氣裝置將一定壓力的氣流輸入到土樣中，且氣壓力逐級(jí)增大．試驗(yàn)過(guò)程中，被驅(qū)趕出來(lái)的水通過(guò)溢流開(kāi)關(guān)流入到量筒中，土樣的體積含水量通過(guò) TDR記錄下來(lái)，壓力表記錄試驗(yàn)過(guò)程中的孔隙氣壓力情況．

2 水-氣二相流數(shù)值模擬模型

?為與多孔介質(zhì)孔隙分布特征有關(guān)的參數(shù)；λ為與曲線形狀有關(guān)的參數(shù)；pmax為最大毛細(xì)壓力．

式(2)用來(lái)表征毛細(xì)壓力-飽和度之間的函數(shù)對(duì)應(yīng)關(guān)系．

VGM模型為

式中：krw和 krg分別為水相和氣相相對(duì)滲透率；τ為迂曲度因子；參數(shù)λ和?之間的函數(shù)關(guān)系滿足λ=1?1?．

式(3)和式(4)分別用來(lái)表征水相和氣相相對(duì)滲透率-飽和度之間的函數(shù)對(duì)應(yīng)關(guān)系．

通過(guò)編制計(jì)算機(jī)程序來(lái)求解上述水-氣二相流模型，由于文章的篇幅有限，詳細(xì)的計(jì)算程序的開(kāi)發(fā)及驗(yàn)證工作見(jiàn)文獻(xiàn)[14-15]．

3 本構(gòu)關(guān)系模型參數(shù)的反演識(shí)別

2.1 基本控制方程

水-氣二相流數(shù)值模擬模型的基本控制方程為水和空氣這2種組分的質(zhì)量守恒方程，即

式中：κ為組分水(w)或空氣(a)；β為水相(l)或氣相(g)；φ為孔隙率；Sβ為β相的飽和度；ρβ為β相流體的密度；Xκβ為κ組分占β相的質(zhì)量分?jǐn)?shù)；q為源匯項(xiàng)；k為固有滲透系數(shù)；rkβ為β相的相對(duì)滲透系數(shù)；βμ為β相的黏滯性系數(shù)；pβ為β相的孔隙壓力；g為重力加速度矢量．

2.2 本構(gòu)關(guān)系模型

在水-氣二相流系統(tǒng)中，VG模型和VGM模型是應(yīng)用最為廣泛、適用性較好的本構(gòu)關(guān)系模型．VG模型為

式中：pcgw為水相和氣相交界面上的毛細(xì)壓力；Swe表示有效水飽和度，Swe= ( Sw? Swr) (Sws?Swr)，Sws為飽和水飽和度，Swr為剩余水飽和度；p0為進(jìn)氣壓力；

水-氣二相流本構(gòu)關(guān)系模型參數(shù)的反演識(shí)別方法：首先進(jìn)行非穩(wěn)態(tài)的水-氣二相流試驗(yàn)，記錄試驗(yàn)過(guò)程中不同時(shí)刻的流場(chǎng)狀態(tài)變量，包括累積出流量、流體飽和度和孔隙壓力等，從而獲取流場(chǎng)狀態(tài)變量的瞬態(tài)測(cè)量值；利用水-氣二相流數(shù)學(xué)模型對(duì)試驗(yàn)進(jìn)行數(shù)值模擬，計(jì)算出與試驗(yàn)過(guò)程相應(yīng)時(shí)刻的流場(chǎng)狀態(tài)變量；采用Marquart-Levenberg最優(yōu)化算法建立最小化目標(biāo)函數(shù)的參數(shù)反演識(shí)別模型，以模擬值與測(cè)量值之間的誤差作為目標(biāo)函數(shù)，通過(guò)模擬值與測(cè)量值的比較，不斷調(diào)整估計(jì)參數(shù)向量，直到滿足收斂條件，從而得到最優(yōu)的水-氣二相流本構(gòu)關(guān)系模型參數(shù)值．本研究所建立的本構(gòu)關(guān)系模型參數(shù)反演識(shí)別流程見(jiàn)圖 2．

本構(gòu)關(guān)系模型參數(shù)的反演識(shí)別是一個(gè)非線性最優(yōu)化問(wèn)題，目標(biāo)函數(shù)為實(shí)測(cè)值與模擬值之間的誤差，加權(quán)最小二乘形式的目標(biāo)函數(shù)為

式中：b為包含 Nc個(gè)參數(shù)的估計(jì)參數(shù)向量，通過(guò)最小化目標(biāo)函數(shù)來(lái)最終確定；oq、oθ和op分別為累積出流量、體積含水量和孔隙氣壓力的觀測(cè)值；qs、θs和ps分別為累積出流量、體積含水量和孔隙氣壓力的模擬值；NQ、NT和 NP分別為累積出流量、體積含水量和孔隙氣壓力的觀測(cè)次數(shù)；wmn(m＝q，θ，p；n=i，j，k)為第n次觀測(cè)時(shí)m變量的權(quán)因子．

圖2 本構(gòu)關(guān)系模型參數(shù)反演識(shí)別流程Fig.2 Flow chart of inverse estimation of parameters under constitutive relationship

將式(7)代入式(6)中，令目標(biāo)函數(shù) ()Ob對(duì)估計(jì)參數(shù)向量b的導(dǎo)數(shù)等于0，可得

在式(6)中，引入迭代指標(biāo) r，用殘差向量(y?ys(b))來(lái)替代y，并且運(yùn)用Marquart-Levenberg優(yōu)化方法[16]，得到迭代形式的估計(jì)參數(shù)向量的增量向量 dr的表達(dá)式為

式中：mr為Marquart參數(shù)；C為對(duì)角比例變換矩陣，=(XTW X)?12；b為N×N階的單位矩陣；?為rCCr阻尼系數(shù)．

在每次迭代過(guò)程中，利用水-氣二相流數(shù)學(xué)模型計(jì)算該迭代步的估計(jì)參數(shù)向量對(duì)應(yīng)的與試驗(yàn)過(guò)程相應(yīng)的模擬值向量；然后將估計(jì)參數(shù)逐個(gè)給定一個(gè)微小的增量，計(jì)算雅克比矩陣和對(duì)應(yīng)的殘差向量；進(jìn)而運(yùn)用Marquart-Levenberg優(yōu)化方法，計(jì)算出估計(jì)參數(shù)向量的增量向量，從而得到下一迭代步的估計(jì)參數(shù)向量；判斷是否收斂，當(dāng)滿足如下收斂條件之一時(shí)迭代停止，即估計(jì)參數(shù)向量中任一參數(shù)變量不大于 2%，或者目標(biāo)函數(shù)的變量不大于 2%，則可得到最優(yōu)的估計(jì)參數(shù)向量．

式(5)可寫成統(tǒng)一的矩陣形式為

式中：e為殘差向量；y為由 ND個(gè)觀測(cè)值組成的向量；W為權(quán)矩陣；s()yb為由 ND個(gè)模擬值組成的向量，可表示為

式中X為 ND×NC階的敏感性(雅克比)矩陣，ND為觀測(cè)次數(shù)，NC為估計(jì)參數(shù)的數(shù)目．組成敏感性矩陣的元素ijX的表達(dá)式為

4 數(shù)值算例

利用水-氣二相流模型對(duì)非穩(wěn)定水-氣二相流試驗(yàn)進(jìn)行數(shù)值模擬，圖3為模擬該試驗(yàn)的計(jì)算模型及網(wǎng)格剖分圖．模型側(cè)面的邊界條件分別為：邊界 1為不透水邊界；上部和下部邊界，即邊界 2和邊界 3，為Dirichlet邊界條件，邊界 2上的氣壓力為大氣壓力，邊界 3上的氣壓力已知，等于壓力控制室內(nèi)的壓力，由壓力表讀出．土柱處于穩(wěn)定飽和狀態(tài)，初始條件由試驗(yàn)給出．

圖3 計(jì)算模型及網(wǎng)格剖分Fig.3 Computational model and meshes

由于試驗(yàn)過(guò)程中，底部氣壓力是逐級(jí)增加的，因此在數(shù)值模擬過(guò)程中要進(jìn)行分段模擬．對(duì)于第 1階段施加的氣壓力，根據(jù)初始已知的初邊值條件進(jìn)行模擬，得到第 1階段末的土柱孔壓、飽和度分布狀態(tài)；開(kāi)始第2階段施加的氣壓力，將第1階段末的滲流狀態(tài)量作為初始條件進(jìn)行模擬，同樣得到第2階段末的滲流狀態(tài)，將其作為模擬第 3階段滲流的初始條件，以此類推，直到模擬完最后一個(gè)階段所施加的氣壓力過(guò)程，計(jì)算出與試驗(yàn)過(guò)程相應(yīng)的模擬值向量．

VG和 VGM 本構(gòu)模型中所包含的參數(shù)有：固有滲透系數(shù) K、孔隙形狀特性參數(shù)λ 和? 、飽和水飽和度 Sws、剩余水飽和度 Swr、迂曲度因子τ、進(jìn)氣壓力p0和最大毛細(xì)壓力 pmax，其中固有滲透系數(shù) K和飽和水飽和度 Sws通過(guò)常水頭試驗(yàn)直接測(cè)定，K＝7.7×10-12,m2，Sws＝0.98，最大毛細(xì)壓力 pmax取試驗(yàn)值，pmax＝1.0×107,N/m2，形狀參數(shù)λ 和? 之間存在函數(shù)關(guān)系為λ=1?1?，因此，反演識(shí)別模型的估計(jì)參數(shù)有4個(gè)：λ、 Swr、τ和 p0．

根據(jù)文獻(xiàn)[17]給出估計(jì)參數(shù)的初始值，利用所建立的反演模型對(duì)估計(jì)參數(shù)進(jìn)行識(shí)別，經(jīng)過(guò) 13次迭代達(dá)到收斂，收斂時(shí)的估計(jì)參數(shù)即為估計(jì)參數(shù)的最優(yōu)值．每次迭代的目標(biāo)函數(shù)值如圖 4所示，估計(jì)參數(shù)的初始值和最優(yōu)值如表1所示；最優(yōu)估計(jì)參數(shù)間的相關(guān)性如表2所示．

表1 估計(jì)參數(shù)的初始值和最優(yōu)值Tab.1 Initial value and optimal value of model estimated parameters

表2 最優(yōu)估計(jì)參數(shù)間的相關(guān)矩陣Tab.2 Correlation matrix of optimal estimated parameters

圖4 每次迭代的目標(biāo)函數(shù)值Fig.4 Objective function value at each iteration

估計(jì)參數(shù)最優(yōu)值對(duì)應(yīng)的水-氣二相流模型的模擬值與試驗(yàn)觀測(cè)值的對(duì)比曲線如圖5～圖7所示，其中圖5為累積出流量的模擬值與實(shí)測(cè)值的對(duì)比，圖6為土柱中壓力表 P3位置(見(jiàn)圖 1)的氣壓力模擬值和實(shí)測(cè)值對(duì)比，圖 7所示的體積含水量為土樣中 T2位置的模擬值和實(shí)測(cè)值．

圖5 累積出流量模擬值與實(shí)測(cè)值的對(duì)比Fig.5 Comparison between simulated and calculated water outflow

圖6 氣壓力模擬值與實(shí)測(cè)值的對(duì)比Fig.6 Comparison between simulated and calculated gas pressure

圖7 體積含水量模擬值與實(shí)測(cè)值的對(duì)比Fig.7 Comparison between simulated and calculated volumetric water content

由此可見(jiàn)，模擬值與實(shí)測(cè)值變化趨勢(shì)基本一致，考慮到試驗(yàn)過(guò)程中一定壓力的空氣入流對(duì)試驗(yàn)土樣的擾動(dòng)性及水-氣二相流模型中假設(shè)土樣為均勻的、各向同性的連續(xù)介質(zhì)所引起的誤差，所得到的擬合結(jié)果能夠較好地說(shuō)明估計(jì)參數(shù)最優(yōu)值的準(zhǔn)確性和參數(shù)反演模型的有效性．

將表 1所示的本構(gòu)關(guān)系模型參數(shù)最優(yōu)值代入到式(4)～式(6)中，即得到與試驗(yàn)中的土樣相對(duì)應(yīng)的毛細(xì)壓力-飽和度、水相和氣相流體的相對(duì)滲透率-飽和度關(guān)系曲線，如圖8和圖9所示，其中圖8中pc單位為 kPa．

圖8 毛細(xì)壓力-飽和度關(guān)系曲線Fig.8 Relationship between capillary pressure and water saturation

圖9 水相和氣相的相對(duì)滲透率-飽和度關(guān)系曲線Fig.9 Relationship between relative permeability and water saturation of water and gas phase

5 結(jié) 語(yǔ)

本研究進(jìn)行了多步的非穩(wěn)定空氣出流試驗(yàn)，測(cè)量了試驗(yàn)過(guò)程中2點(diǎn)的孔隙壓力、氣相飽和度及累積出流量的變化過(guò)程；建立了考慮水相、氣相流動(dòng)及相互溶混的水-氣二相流數(shù)值模擬模型，采用 Marquart-Levenberg最優(yōu)化算法建立了水-氣二相流本構(gòu)關(guān)系模型參數(shù)反演識(shí)別模型，結(jié)合非穩(wěn)定水-氣二相流試驗(yàn)，對(duì)水-氣二相流的本構(gòu)關(guān)系模型(VG 和 VGM)參數(shù)進(jìn)行反演識(shí)別．得到水-氣二相流本構(gòu)關(guān)系模型參數(shù)的最優(yōu)估計(jì)值，分析了目標(biāo)函數(shù)值在迭代過(guò)程中的變化過(guò)程，統(tǒng)計(jì)了最優(yōu)估計(jì)參數(shù)間的相關(guān)系數(shù)矩陣，對(duì)所建立的參數(shù)反演識(shí)別模型進(jìn)行了評(píng)價(jià)，將參數(shù)最優(yōu)估計(jì)值對(duì)應(yīng)的模擬值與試驗(yàn)過(guò)程的實(shí)測(cè)值進(jìn)行對(duì)比分析，二者吻合較好，證明了估計(jì)參數(shù)最優(yōu)值的準(zhǔn)確性和參數(shù)反演識(shí)別模型的有效性．

在地下水多相流數(shù)值模擬研究中，確定水-氣二相流體的本構(gòu)關(guān)系模型參數(shù)是最基本、也是最關(guān)鍵的．結(jié)合非穩(wěn)定水-氣二相流試驗(yàn)，通過(guò)建立參數(shù)反演識(shí)別模型來(lái)確定水-氣二相流本構(gòu)關(guān)系的方法是很有優(yōu)越性的，盡管該方法在具體實(shí)施過(guò)程中還可能會(huì)遇到許多困難，但是隨著試驗(yàn)技術(shù)水平的提高及模擬試驗(yàn)的數(shù)學(xué)模型的進(jìn)一步完善，所確定的水-氣二相流模型的精度將會(huì)繼續(xù)提高，因此本文的研究成果可以為確定多相流系統(tǒng)中與各相流體(包括水相、氣相和非水相流體（NAPL)）有關(guān)的本構(gòu)關(guān)系提供幫助和技術(shù)支持．

［1］ 林學(xué)鈺. “地下水科學(xué)與工程”學(xué)科形成的歷史沿革及其發(fā)展前景[J]. 吉林大學(xué)學(xué)報(bào)：自然科學(xué)版，2007，37(2)：209-215.Lin Xueyu. Historical change and prospect of discipline development of "Groundwater Science and Engineering"[J]. Journal of Jilin University：Earth Science Edition，2007，37(2)：209-215(in Chinese).

［2］ Helmig R. Multiphase Flow and Transport Processes in the Subsurface，a Contribution to the Modeling of Hydrosystems [M]. Berlin，Heidelberg：Springer-Verlag，1997.

［3］ Brooks R H，Corey A T. Hydraulic properties of porous medium [J]. Hydrology Paper of Colorado State University Fort Collins，1964(3)：1-27.

［4］ van Genuchten M T. A closed-form equation for predict-ing the hydraulic conductivity of unsaturated soils [J].Soil Science Society America Journal，1980，44：892-898.

［5］ Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media [J]. Water Resources Research，1976，12(3)：513-522.

［6］ Leverett M C. Capillary behavior in porous solids [J].Trans AIME，1941，142：152-169.

［7］ Parker J C，Lenhard R J，Kuppusamg T. A parametric model for constitutive properties governing multi-phase flow in porous media [J]. Water Resources Research，1987，23：618-624.

［8］ Dane J H，Oostrom M，Missildine B C. An improved method for the determination of capillary pressuresaturation curves involving TCE，water and air [J].Journal of Contaminant Hydrology，1992，11(1/2)：69-81.

［9］ Chen J，Hopmans J W，Grismer M E. Parameter estimation of two-fluid capillary pressure-saturation and permeability functions [J]. Advances in Water Resources，1999，22(5)：479-493.

［10］ 薛 強(qiáng)，馮夏庭，梁 冰，等. 水氣兩相流系統(tǒng) K-SP模型參數(shù)反演的最優(yōu)估計(jì)[J]. 水科學(xué)進(jìn)展，2005，16(4)：488-495.Xue Qiang，F(xiàn)eng Xiating，Liang Bing，et al. Optimal estimation of parameters inversion for the relationship of permeability-saturation-pressure in water-air phase flow system [J]. Advances in Water Science，2005，16(4)：488-495(in Chinese).

［11］ Gardner W R. Calculation of capillary conductivity from pressure plate outflow data [J]. Soil Science Society America Journal，1956，20：317-320.

［12］ 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 America Journal，1994，58：647-652.

［13］ Eching S O，Hopmans J W. Optimization of hydraulic functions from transient outflow and soil water pressure data [J]. Soil Science Society America Journal，1993，57：1167-1175.

［14］ 孫冬梅. 水-氣二相流飽和-非飽和滲流場(chǎng)分析及其應(yīng)用研究[D]. 南京：河海大學(xué)水利水電工程學(xué)院，2007.Sun Dongmei. Study on Saturated-Unsaturated Seepage Based on Water-Air Two-Phase Flow Theory and Applications [D]. Nanjing：College of Water Conservancy and Hydropower Engineering，Hohai University，2007(in Chinese).

［15］ 孫冬梅，朱岳明，張明進(jìn). 非飽和帶水-氣二相流數(shù)值模擬研究[J]. 巖土工程學(xué)報(bào)，2007，29(4)：560-565.Sun Dongmei，Zhu Yueming，Zhang Mingjin. Study on numerical model for water-air two-phase flow in unsaturated soil [J]. Chinese Journal of Geotechnical Engineering，2007，29(4)：560-565(in Chinese).

［16］ Poeter E P，Hill M C. Inverse models：A necessary next step in groundwater modeling [J]. Ground Water，1997，35(2)：250-260.

［17］ Arya L M，Paris J F. A physicoempirical model to predict the soil moisture characteristic from particle-size distribution and bulk density data [J]. Soil Science Society America Journal，1981，45(6)：1023-1030.