滯后現(xiàn)象對(duì)CO2深部咸水層地質(zhì)封存的影響

郭朝斌，張可霓，魯維豐，凌璐璐，郭振東

（1.北京師范大學(xué)水科學(xué)研究院，北京 100875；2.神華鄂爾多斯煤制油分公司，內(nèi)蒙古 鄂爾多斯伊旗 017209）

溫室氣體排放增加導(dǎo)致溫室效應(yīng)加劇，目前全球變暖已經(jīng)成為國(guó)際社會(huì)關(guān)注的熱點(diǎn)問(wèn)題。CO2的封存問(wèn)題就是在這樣的背景下提出來(lái)的，CO2的捕獲和封存技術(shù)目前已經(jīng)成為一種國(guó)際公認(rèn)的減排措施，被認(rèn)為有巨大的減排潛力［1］。地下咸水層在世界各地廣泛分布，而且規(guī)模大，是切實(shí)可行和最具發(fā)展前景的減排途徑，被認(rèn)為是長(zhǎng)期地質(zhì)封存的最有效方法之一［2］。

近年來(lái)國(guó)內(nèi)外很多學(xué)者對(duì)CO2地質(zhì)封存的數(shù)值模擬陸續(xù)展開(kāi)了研究。R.Juanes等［3］通過(guò)對(duì)CO2捕獲機(jī)制和相對(duì)滲透率的研究及數(shù)值模擬，表明交替注入水和CO2可以提高注入效率。注入的水使大面積CO2暈分裂，增大吸濕過(guò)程的時(shí)間和面積從而增強(qiáng)CO2的捕獲和固定，但這種方法會(huì)導(dǎo)致井孔底部壓力增大，受到蓋層的封閉性能、管理和經(jīng)濟(jì)方面等的限制。理解大體積CO2暈的流體力學(xué)行為對(duì)提高咸水層中CO2地質(zhì)封存技術(shù)非常重要，Christine Doughty［4］對(duì)美國(guó)加州圣華金河峽谷南部Kimberlina發(fā)電廠進(jìn)行數(shù)學(xué)建模研究，利用TOUGH2對(duì)傾斜地層中CO2地質(zhì)封存大規(guī)模初步試驗(yàn)中CO2暈的演變進(jìn)行模擬，結(jié)果表明在4年注入周期后，注入的100×104t CO2中約20%為溶液相，25年后，CO2暈達(dá)到穩(wěn)定，溶液相態(tài)比例上升到38%，62%為超臨界狀態(tài)，其中只有3%為自由態(tài)。隨著時(shí)間繼續(xù)，暈的演變僅表現(xiàn)為緩慢的溶解作用，并伴有飽和溶解CO2咸水的對(duì)流。國(guó)內(nèi)對(duì)提高CO2封存能力、安全評(píng)估等方面研究中，許雅琴等［5］通過(guò)模擬提高低滲地層中滲透率等探討如何提高咸水層CO2封存注入率，但對(duì)滯后效應(yīng)在CO2封存方面的研究相對(duì)較少。

1 方法和理論

1.1 超臨界狀態(tài)CO2的吸濕和疏干過(guò)程

將超臨界狀態(tài)CO2封存到深部咸水層會(huì)涉及到諸多多相流問(wèn)題，注入的CO2為非潤(rùn)濕相，水及少量溶解到水中的CO2作為潤(rùn)濕相。由于潤(rùn)濕相和非潤(rùn)濕相流體性質(zhì)不同，運(yùn)移過(guò)程中因接觸角不同和殘余氣相的產(chǎn)生會(huì)對(duì)流體運(yùn)移產(chǎn)生影響。注入時(shí)CO2暈連續(xù)增長(zhǎng)，超臨界狀態(tài)CO2(表征為氣相)作為非潤(rùn)濕相驅(qū)替含水層中潤(rùn)濕相(咸水)，即疏干過(guò)程。對(duì)于隨后停止注入階段，當(dāng)CO2暈在浮力作用下向上運(yùn)移(圖1)，不同位置經(jīng)歷過(guò)程不同，在運(yùn)移前端為疏干過(guò)程，在暈尾端，潤(rùn)濕相驅(qū)替非潤(rùn)濕相，為吸濕過(guò)程，有部分殘余CO2被捕獲，變?yōu)椴豢蛇\(yùn)動(dòng)的殘余氣相，進(jìn)而影響相對(duì)滲透率和毛細(xì)壓力的計(jì)算。不同位置和過(guò)程中殘余氣相不同，使得相同飽和度下相對(duì)滲透率和毛細(xì)壓力不相同，在特征曲線上表現(xiàn)出滯后現(xiàn)象。

圖1 CO2向上運(yùn)移過(guò)程暈后邊緣殘余示意圖［3］Fig.1 Schematic diagram showing the trail of residual CO2as the plume migrates upward

1.2 數(shù)值模擬方法

數(shù)學(xué)模型和數(shù)值模擬在評(píng)估地質(zhì)封存CO2可行性上具有重要的作用，它們是設(shè)計(jì)和實(shí)施CO2地質(zhì)處置必須的工具［6］?？梢栽u(píng)估預(yù)測(cè)地層的封存能力，評(píng)估CO2地質(zhì)儲(chǔ)存的可行性和可靠性，解釋CO2在地層中的運(yùn)移變化行為以及在注入到地層中后對(duì)其運(yùn)移過(guò)程進(jìn)行監(jiān)測(cè)和分析。

在大多數(shù)通用多相流數(shù)值模擬軟件中，相對(duì)滲透率和毛細(xì)壓力的計(jì)算函數(shù)僅與當(dāng)前飽和度有關(guān)，然而實(shí)驗(yàn)表明［7］，相對(duì)滲透率和毛細(xì)壓力不僅和當(dāng)前區(qū)域的飽和度有關(guān)，還和流體經(jīng)歷的過(guò)程及歷史飽和度有關(guān)。盡管滯后效應(yīng)特征曲線在石油工業(yè)已經(jīng)被應(yīng)用，但是在包氣帶、地?zé)醿?chǔ)藏工程、核廢料地質(zhì)處理以及CO2地質(zhì)封存等這些兩相流的問(wèn)題中很少被應(yīng)用到［8～9］。

目前在CO2地質(zhì)封存方面，應(yīng)用較為普遍的模擬軟件為T(mén)OUGH2及2008年正式公開(kāi)發(fā)布用于大規(guī)模并行計(jì)算的并行版本TOUGH2-MP［10］。目前成功應(yīng)用最大的模型是東京灣的CO2地質(zhì)處置模型，其單元數(shù)超過(guò)1 千萬(wàn)個(gè)［11］。TOUGH+ 是繼 TOUGH2［12］之后的新一代TOUGH家族的模擬器，在繼承了TOUGH2所有功能外，主要在數(shù)組處理效率方面提高模擬的計(jì)算效率，同時(shí)也提高了一些模擬功能。TOUGH+CO2是TOUGH+系列中模擬CO2封存的模塊［13］，最新版本TOUGH+CO2中，在 iTOUGH2［14］中關(guān)于滯后效應(yīng)模塊的基礎(chǔ)上改進(jìn)更新了滯后效應(yīng)模擬的功能。

1.3 TOUGH+CO2中相對(duì)滲透率相關(guān)理論

滯后效應(yīng)特征曲線不僅取決于當(dāng)前飽和度，還和歷史飽和度以及經(jīng)歷的過(guò)程即疏干過(guò)程(非潤(rùn)濕相替代潤(rùn)濕相)或吸濕過(guò)程(潤(rùn)濕相替換非潤(rùn)濕相)相關(guān)［15～17］。

CO2深部咸水層封存中主要涉及多相流系統(tǒng)為H2O－CO2－NaCl系統(tǒng)，超臨界狀態(tài)CO2在TOUGH+CO2中表征為非水相，即氣相，溶解于水中的CO2表征為溶液相；計(jì)算系統(tǒng)中流體相的相對(duì)滲透率的常用方法為 van Genuchten［18］公式。

考慮相對(duì)滲透率滯后現(xiàn)象時(shí)，TOUGH+CO2采用Parker和Lenhard根據(jù)無(wú)滯后效應(yīng)van Genuchten方程修改得到的公式［14］，

式中:krl——液相相對(duì)滲透率，無(wú)量綱；

krg——?dú)庀嘞鄬?duì)滲透率，無(wú)量綱；

式中:Sl——液相飽和度，無(wú)量綱；

Slr——液相殘余飽和度，無(wú)量綱；

Sls——飽和水飽和度，無(wú)量綱；

Sgr——?dú)堄鄽庀囡柡投?，無(wú)量綱。

圖2 歷史飽和度對(duì)相對(duì)滲透率的影響Fig.2 Effect of the history saturation on relative permeability

1.4 TOUGH+CO2中毛細(xì)壓力相關(guān)理論

在計(jì)算H2O－CO2系統(tǒng)中流體相的毛細(xì)壓力時(shí)經(jīng)常用到的方法為 van Genuchten［18］函數(shù)。

式中:Pcap——毛細(xì)壓力(Pa)；

P0——進(jìn)入毛細(xì)壓力(Pa)。

考慮毛細(xì)壓力滯后現(xiàn)象時(shí)采用基于van Genuchten方程修改得到的公式［14］，

式中:p——流體流動(dòng)過(guò)程，疏干過(guò)程(d)或吸濕過(guò)程(w)；

mp——函數(shù)適配參數(shù) mp=(nγ－1)/nγ，n、γ 為適配參數(shù)。

2 應(yīng)用

2.1 模型基本參數(shù)

通過(guò)對(duì)案例模型的模擬，可以分析相對(duì)滲透率和毛細(xì)壓力滯后效應(yīng)對(duì)CO2咸水層地質(zhì)封存的影響，為了研究方便，假設(shè)模型為理想條件模型。

進(jìn)行CO2深部咸水層封存模擬時(shí)，需要給出的初始條件有壓力、溫度、鹽度以及流體中CO2的質(zhì)量分?jǐn)?shù)。模型模擬的咸水層在地下1200m深，模擬范圍為1km×1km，模型厚度為100m，模型上下邊界為低滲透地層，其中蓋層為5m厚泥巖，注入層為85m厚砂巖，下伏層為10m厚泥巖(圖3)。模型共剖分成64000個(gè)網(wǎng)格，其中注入井附近網(wǎng)格進(jìn)行加密處理。四周邊界條件視為無(wú)流量邊界，對(duì)于封閉系統(tǒng)，即儲(chǔ)層被完全封閉，其CO2儲(chǔ)存量由在儲(chǔ)層中壓力積聚引起孔隙體積的膨脹和咸水密度的增加來(lái)得到。

圖3 模型范圍XZ剖面示意圖Fig.3 Profile of the model in the XZ direction

此模型中只考慮等溫過(guò)程，初始溫度為37℃，模型的頂部邊界壓力為1.20×107Pa，初始?jí)毫Ψ植几鶕?jù)靜水壓力平衡計(jì)算得到。超臨界狀態(tài)下CO2的性質(zhì)根據(jù)實(shí)驗(yàn)數(shù)據(jù)加入到程序，其他關(guān)鍵物理參數(shù)見(jiàn)表1。

表1 模型采用關(guān)鍵物理參數(shù)Table 1 Key properties of rocks in the model

方案1:不考慮滯后現(xiàn)象對(duì)CO2封存的影響，在計(jì)算中，相對(duì)滲透率和毛細(xì)壓力的計(jì)算方式均選擇van Genuchten函數(shù)。模型區(qū)域中介質(zhì)的殘余氣相飽和度Sgr為0.177，整個(gè)含水層中的殘余氣相飽和度均為固定Sgr。在數(shù)值模擬中，根據(jù)固定值Sgr，計(jì)算相對(duì)滲透率特征曲線，在驅(qū)替過(guò)程中(包括疏干和吸濕過(guò)程)，相對(duì)滲透率均沿此曲線變化，同一流體飽和度下對(duì)應(yīng)相同相對(duì)滲透率。毛細(xì)壓力的計(jì)算亦是如此。

方案2:考慮滯后現(xiàn)象對(duì)CO2封存的影響。相對(duì)滲透率和毛細(xì)壓力計(jì)算函數(shù)選擇修改后的公式。模型區(qū)域中殘余氣相飽和度的最大可能值Sgrmax為0.177。在數(shù)值計(jì)算中，根據(jù)驅(qū)替過(guò)程轉(zhuǎn)換點(diǎn)的氣相飽和度和經(jīng)歷的過(guò)程計(jì)算相應(yīng)的殘余氣相飽和度，然后通過(guò)插值得到掃描曲線，與主線不同。殘余液相飽和度Slr為0.2，krgmax為 1。

CO2注入以0.926 kg/s(約3×104t/a)恒定速率注入4年，為觀察停止注入后的運(yùn)移狀態(tài)，在第5～50年停止注入，模擬時(shí)間總共為50年。

2.2 結(jié)果分析

2.2.1 對(duì)封存量的影響

在模擬時(shí)間為50年時(shí)，99%的CO2存儲(chǔ)于注入層中，只有1%存在蓋層和底部地層中(表2)。因浮力作用CO2向上運(yùn)移，由于蓋層滲透率低，封閉性好，CO2大部分被封存在注入層中。所以兩個(gè)方案中封存的CO2總量基本一致，但CO2存在相態(tài)比例不同。

表2 模擬時(shí)間為50年時(shí)各層CO2含量Table 2 Mass of CO2in various forms at 50 years (kg)

2.2.2 注入CO2的存在形式及其分布

圖4～5為兩種方案中封存CO2相態(tài)比例變化。氣相是指整個(gè)地層中超臨界狀態(tài)CO2，包含可移動(dòng)的自由態(tài)CO2和捕獲封存不可移動(dòng)氣相CO2。隨著模擬時(shí)間的增加，CO2和咸水層接觸時(shí)間增加，氣相CO2逐漸溶解于溶液中成為溶液相。

圖4 無(wú)滯后效應(yīng)模型整個(gè)地層中各相態(tài)CO2質(zhì)量分?jǐn)?shù)Fig.4 Mass fraction of CO2in the non-hysteretic case over the entire model

圖5 滯后效應(yīng)模型整個(gè)地層中各相態(tài)CO2質(zhì)量分?jǐn)?shù)Fig.5 Mass fraction of CO2in the hysteretic case over the entire model

模擬到50年時(shí)，無(wú)滯后方案中氣相形態(tài)逐漸減少到70%，滯后方案中氣相保持約45%的比例。產(chǎn)生這種差異的原因是無(wú)滯后方案中，殘余氣相飽和度為0.177，即含水層介質(zhì)都默認(rèn)捕獲此數(shù)值殘余氣相飽和度，只有氣相飽和度高于此殘余氣相飽和度時(shí)才運(yùn)移。而在滯后方案中，是可能存在的最大氣相飽和度，即每個(gè)網(wǎng)格元素中根據(jù)經(jīng)歷的過(guò)程不同，轉(zhuǎn)換點(diǎn)飽和度不同，殘余氣相飽和度不同。因此滯后效應(yīng)方案中捕獲的殘余氣相總是小于等于無(wú)滯后效應(yīng)方案中的殘余氣相，因此表現(xiàn)出在模擬時(shí)間到50年時(shí)，無(wú)滯后效應(yīng)方案中封存的氣相CO2比例大于滯后效應(yīng)中的氣相比例。同樣，封存的CO2在含水層中的分布也不同，圖6為在CO2注入階段，即在疏干過(guò)程中氣相CO2飽和度的分布。在疏干過(guò)程中，兩種方案的相對(duì)滲透率特征曲線均沿主線變化，不涉及滯后現(xiàn)象，但因?yàn)闅堄鄽庀囡柡投鹊牟煌鄳?yīng)特征曲線則不同。以注入點(diǎn)附近區(qū)域?yàn)槔瑹o(wú)滯后效應(yīng)的氣相krg總小于滯后效應(yīng)方案的krg，因此在暈圖上表現(xiàn)出無(wú)滯后效應(yīng)方案中氣相CO2運(yùn)移緩慢(圖7)。

第4年末CO2停止注入后，在暈的尾端發(fā)生吸濕過(guò)程。如圖7，無(wú)滯后方案捕獲更多的殘余氣相，故CO2暈分布集中在運(yùn)移過(guò)程中，在滯后方案中部分氣相CO2則快速運(yùn)移到蓋層下，由于蓋層的封閉性好，氣相CO2暈逐漸向周?chē)鷶U(kuò)散。

圖8表示數(shù)值計(jì)算中相對(duì)滲透率特征曲線。在無(wú)滯后方案中，特征曲線沿著疏干過(guò)程逆向變化，即相互驅(qū)替過(guò)程只有一條曲線。而在滯后方案中，相對(duì)滲透率和毛細(xì)壓力不僅與當(dāng)前區(qū)域的飽和度有關(guān)，還和經(jīng)歷的過(guò)程和歷史飽和度有關(guān)。例如，在注入點(diǎn)網(wǎng)格區(qū)域，疏干過(guò)程中Sg增加，相應(yīng)的krg增加，krl減少到0。在隨后停止注入的吸濕過(guò)程中，Sg逐漸減少，krg減少，krl逐漸增大。但在同一Sg時(shí)，兩個(gè)過(guò)程的相對(duì)滲透率不同，產(chǎn)生滯后的現(xiàn)象。同樣，在中間某網(wǎng)格處，相對(duì)滲透率變化趨勢(shì)相同，只是轉(zhuǎn)換點(diǎn)不同，進(jìn)而相對(duì)滲透率變化曲線不同。

圖6 CO2注入過(guò)程中無(wú)滯后效應(yīng)方案(上)和滯后效應(yīng)方案(下)中氣相飽和度分布圖Fig.6 CO2plume evolution during the injection period in the non-hysteretic case(top)and the hysteretic case(bottom)

圖7 停止注入后無(wú)滯后效應(yīng)方案(上)和滯后效應(yīng)(下)方案中氣相飽和度隨時(shí)間變化分布圖Fig.7 CO2plume evolution during the post-injection period in the non-hysteretic case(top)and the hysteretic case(bottom)

圖8 模型中注入點(diǎn)網(wǎng)格相對(duì)滲透率變化曲線Fig.8 Relative permeability path for the injection grid

圖9～10對(duì)應(yīng)兩種方案中毛細(xì)壓力特征曲線，選取模型中注入點(diǎn)網(wǎng)格和中間網(wǎng)格兩個(gè)不同位置得到其毛細(xì)特征曲線。滯后方案為一次完整疏干和吸濕過(guò)程，在中間網(wǎng)格處，氣相飽和度從0逐漸增大到轉(zhuǎn)換點(diǎn)飽和度，然后沿掃描曲線減少到相應(yīng)殘余氣相飽和度，并未沿疏干過(guò)程主線變化，導(dǎo)致在同一氣相飽和度下毛細(xì)壓力不同。無(wú)滯后模型中，毛細(xì)壓力均沿主線變化。在實(shí)際情況中，僅用實(shí)驗(yàn)測(cè)試得到的殘余氣相飽和度賦值給含水層來(lái)進(jìn)行數(shù)值模擬，在封存能力、安全評(píng)估等方面會(huì)造成誤差。

圖9 滯后效應(yīng)模型中不同位置毛細(xì)壓力變化曲線Fig.9 Capillary pressure path for grids in the hysteretic

圖10 無(wú)滯后效應(yīng)模型不同位置毛細(xì)壓力變化曲線Fig.10 Capillary pressure path for grids in the non-hysteretic case

3 結(jié)論

在CO2咸水層封存的數(shù)值模擬計(jì)算中，相對(duì)滲透率和毛細(xì)壓力計(jì)算函數(shù)是重要計(jì)算參數(shù)之一。通過(guò)設(shè)計(jì)對(duì)應(yīng)無(wú)滯后效應(yīng)和考慮相對(duì)滲透率、毛細(xì)壓力滯后的兩種方案，運(yùn)用數(shù)值模擬軟件TOUGH+CO2進(jìn)行模擬，結(jié)果表明不同相對(duì)滲透率和毛細(xì)壓力計(jì)算函數(shù)的選擇表現(xiàn)出較大差異。滯后效應(yīng)方案能更詳細(xì)描述每個(gè)網(wǎng)格中的相對(duì)滲透率和毛細(xì)壓力變化情況。在蓋層封閉性較好時(shí)，滯后效應(yīng)對(duì)CO2封存總量沒(méi)有影響，但相態(tài)比例不同。

不同相對(duì)滲透率和毛細(xì)壓力計(jì)算函數(shù)的選擇表現(xiàn)出CO2暈運(yùn)移分布不同，對(duì)實(shí)際封存方案設(shè)計(jì)，封存能力評(píng)估以及環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)等都會(huì)帶來(lái)不同影響，因此在CO2咸水層封存數(shù)值模擬中，尤其是停止注入之后應(yīng)引入與滯后現(xiàn)象相關(guān)的模擬計(jì)算。

［1］郝艷軍，楊頂輝.二氧化碳地質(zhì)封存問(wèn)題和地震監(jiān)測(cè)研究進(jìn)展［J］.地球物理學(xué)進(jìn)展，2012，27(6):2369 － 2383.［HAO Y J，YANG D H.Research progress of carbon dioxide capture and geological sequestration problem and seismic monitoring research［J］.Progress in Geophys，2012，27(6):2369－2383.(in Chinese)］

［2］張煒，呂鵬.二氧化碳地質(zhì)封存中對(duì)流混合過(guò)程的研究進(jìn)展［J］.水文地質(zhì)工程地質(zhì)，2013，40(2):101－107.［ZHANG W，LV P.Density-driven convection in carbon dioxide geological storage:a review ［J］.Hydrogeology ＆ Engineering Geology，2013，40(2):101－107.(in Chinese)］

［3］Juanes R，Spiteri E J，Orr F M，et al.Impact of relative permeability hysteresis on geological CO2storage［J］.Water Resource Research，2006，42(12):1－13.

［4］Doughty C.Investigation of CO2plume behavior for a large-scale pilot test of geologic carbon storage in a saline formation ［J］.Transport in porous media，2010，82(1):49 －76.

［5］許雅琴，張可霓，王洋.基于數(shù)值模擬探討提高咸水層CO2封存注入率的途徑［J］.巖土力學(xué)，2012，33(12):3825－3832.［XU Y Q，ZHANG K N，WANG Y.Numerical investigation for enhancing injectivity of CO2storage in saline aquifers［J］.Rock and Soil Mechanics，2012，33(12):3825 －3832.(in Chinese)］

［6］張曉宇，成建梅，劉軍，等.CO2地質(zhì)處置研究進(jìn)展［J］.水文地質(zhì)工程地質(zhì)，2006，33(4):85－89.［ZHANG X Y，CHENG J M，LIU J，et al.An overview of underground sequestration of carbon dioxide［J］.Hydrogeology ＆ Engineering Geology，2006，33(4):85 －89.(in Chinese)］

［7］周琦，朱學(xué)謙，劉傳喜.毛細(xì)管力滯后曲線和相滲透率滯后曲線及處理［J］.油氣采收率技術(shù)，1999，6(3):47－50.［ZHOU Q，ZHU X Q，LIU C X.Curves of hysteretic phenomenon in relative permeability and capillary pressure［J］.Oil＆ Gas Recovery Technology，1999，6(3):47－50.(in Chinese)］

［8］Finsterle S.Multiphase inverse modeling:Review and iTOUGH2 applications ［J］.Vadose Zone Journal，2004，3(3):747－762.

［9］Zhang Y Q，Oldenburg C M，F(xiàn)insterle S，et al.Systemlevel modeling for economic evaluation of geological CO2storage in gas reservoirs［J］.Energy Conversion and Management，2007，48(6):1827 －1833.

［10］Zhang K，Doughty C，Wu Y S，et al.Efficient parallel simulation of CO2geologic sequestration in saline aquifers［C］//Lawrence Berkeley National Laboratory.SPE Reservoir Simulation Symposium.Berkeley:Society of Petroleum Engineers，2007.

［11］Yamamotoa H，Zhang K，Karasakib K，et al.Largescale numerical simulation of CO2geologic storage and its impact on regional groundwater flow:A hypothetical case study at Tokyo Bay［J］.Japan Energy Procedia，2009，1(1):1871－1878.

［12］Pruess K，Oldenburg C，Moridis G.TOUGH2 User’s Guide［R］. LBNL-43134.Berkeley:Lawrence Berkeley National Laboratory，1999.

［13］Zhang K，Moridis G，Pruess K.TOUGH+CO2:A multiphase fluid-flow simulator for CO2geologic sequestration in saline aquifers［J］.Computers ＆Geosciences，2011，37(6):714 －723.

［14］Finsterle S.iTOUGH2 user’s guide［R］.LBNL－40040.Berkeley:Lawrence Berkeley National Laboratory，1999.

［15］Doughty C. Modeling geologic storage of carbon dioxide Comparison of non-hysteretic and hysteretic characteristic curves［J］.Energy Conversion and Management，2007，48(6):1768 －1781.

［16］王大純，張人權(quán)，史毅虹，等.水文地質(zhì)學(xué)基礎(chǔ)［M］.北京:地質(zhì)出版社，1995.

［17］薛禹群.地下水動(dòng)力學(xué)［M］.北京:地質(zhì)出版社，1997.

［18］Van Genuchten M Th，A closed-form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Science Society of America Journal，1980，44(5):892－898.