佘冬立 韓 笑 孫梟沁 唐勝強 王洪德

(河海大學農(nóng)業(yè)科學與工程學院, 南京 210098)

0 引言

土壤孔隙結構是土壤中水、氣、養(yǎng)分等傳輸和貯存的場所[1]。在不考慮固液界面上其它復雜的物理化學過程時,土壤水分入滲等水力性能主要受到孔隙結構的影響[2-3]?？紫段⒂^結構通常能很好地解釋宏觀性能,孔隙尺度上的滲流模擬也成為聯(lián)系微觀結構和宏觀性能的有效途徑[4-5]。精準表征孔隙內(nèi)部微觀結構,是預測土壤水分滲流等過程的基礎。近年來,計算機斷層掃描(CT)技術由于可以無擾動、多方位地獲取土壤結構的無損圖像,成為精細化獲取土壤孔隙結構的可靠手段[6-9]。利用CT掃描技術,結合具有孔隙拓撲結構的網(wǎng)絡模型,能夠在孔隙尺度上模擬和預測多孔介質(zhì)的滲流能力[2,10]。劉向君等[11]基于微CT技術,結合圖像處理等算法,構建了等效孔隙網(wǎng)絡模型,實現(xiàn)了孔隙尺度的滲流模擬并計算得到砂巖的絕對滲透率。閆國亮等[12]則利用過程模擬法和改進的最大球算法,結合逾滲理論開展儲層巖石的孔隙結構研究及其滲透率模擬。FAN等[13]利用微CT掃描和Avizo軟件,建立了煤的等效孔隙網(wǎng)絡模型,并用Comsol軟件模擬了單相水滲流過程。WEI等[14]利用等效孔隙網(wǎng)絡模型定量分析了馬蘭黃土3種類型土壤的內(nèi)部微觀結構,并討論了微觀結構對土壤滲透率的影響。LI等[15]同樣以西北馬蘭黃土為研究對象,著重模擬研究了土壤滲透性能的各向異性。但現(xiàn)階段該方法在土壤方面的應用多集中在黃土高原地區(qū),在濱海土壤方面的應用還較少,而濱海圍墾區(qū)土壤會受到海水浸漬的影響,土壤和地下水接觸過程中會發(fā)生一系列的物理化學反應,使土壤的空間結構呈現(xiàn)出特殊性[16]。

本文基于CT掃描技術,結合孔隙網(wǎng)絡模型,對沿海地區(qū)不同圍墾區(qū)剖面土壤的孔隙結構進行統(tǒng)計分析,并模擬單相水的滲流過程,以利于從微觀角度分析孔隙分布及拓撲關系,深入了解沿海圍墾區(qū)土壤中水分滲流特征及其主要影響因素。

1 材料與方法

1.1 試驗區(qū)概況及土樣采集

試驗于2019年9—10月開展。試驗樣品取自江蘇省沿海地區(qū),包括如東、濱海和東臺3個縣市(圖1,圖中Y1～Y3表示3個時期的圍墾區(qū),Y1為2000年之后,Y2為1960—2000年,Y3為1940—1960年)。其中東臺市試驗區(qū)位于東臺市東北角(32°33′～32°57′N,120°53′～121°07′E),區(qū)域年平均氣溫14～15℃,年平均降水量1 059.8 mm,降水量季節(jié)變化明顯;如東縣試驗區(qū)位于如東縣東北部(32°30′～32°40′N,121°00′～120°11′E),瀕臨黃海,年平均降水量為1 026 mm,平均氣溫15℃;濱?？h試驗區(qū)位于鹽城市東北部(33°43′～34°23′N,119°37′～120°20′E),東瀕黃海,平均氣溫為14.4℃,年平均降水量964.8 mm。

圖1 取樣點位置圖

各墾區(qū)均選擇10 m×10 m樣地作為每個研究區(qū)域的采樣點,樣地均采用稻麥輪作模式,采樣時水稻處于黃熟期。利用PVC柱(直徑5.0 cm,高度5.0 cm)采集樣地0～20 cm、20～40 cm和40～60 cm土層原狀土樣,采用CT掃描獲得孔隙結構圖像。PVC柱在運輸過程中用塑料薄膜包裹,置于防震包裝袋中,以避免土壤擾動。采用標準環(huán)刀(直徑5.0 cm,高度5.0 cm)采集各土層原狀土樣,利用定水頭法測定土壤飽和導水率Ks。

1.2 CT掃描及圖像處理

采用工業(yè)CT掃描儀(Phenix Nano tom S,GE,美國)進行圖像掃描,像素分辨率為25 μm,設置參數(shù)包括掃描電壓100 kV、掃描電流100 μA和曝光時間1.25 s。CT掃描各樣本柱獲約2 018幅二維切片圖像。選用中值濾波算法對圖像進行降噪處理。圖像分割是CT技術應用中獲取土壤孔隙結構的重要步驟[6],在本研究中,通過對比處理后的圖像與原始圖像,根據(jù)實際孔隙條件手動確定閾值,獲得最終的二值圖像。

1.3 表征單元體確定

為保證提取的結構網(wǎng)絡具有代表性,同時提高建模和計算效率,有必要確定表征單元體(REV)[17]。本文選取孔隙度作為確定REV尺寸的重要參數(shù)和約束條件[5]。當孔隙度變化受立方體尺寸影響不明顯時,即可確定為REV的合理尺寸。通過對27組土樣的REV進行“孔隙度-立方體尺寸”分析,得到REV尺寸為300體素×300體素×300體素。

1.4 孔隙網(wǎng)絡模型

采用最大球算法提取孔隙網(wǎng)絡模型[18]。孔隙網(wǎng)絡模型將土壤孔隙結構簡化為孔隙和喉道兩部分,通過對網(wǎng)絡結構參數(shù)定量分析,確定孔隙網(wǎng)絡連通性,最終實現(xiàn)土壤滲透能力的模擬。其中不連通的“死孔”在模擬過程中不起作用,為提高模擬效率,在提取網(wǎng)絡模型時,只保留土壤中連通孔隙。通過對土樣REV分別提取孔隙網(wǎng)絡,土壤孔隙結構被簡化成一個幾何結構(圖2),其中球代表孔隙,棍則代表孔隙之間的喉道。

圖2 采集土樣提取的孔隙網(wǎng)絡模型

1.5 單相水滲流模擬

絕對滲透率用于衡量飽和單相流體通過土壤孔隙的能力[11]。利用構建的孔隙網(wǎng)絡結構,選用常態(tài)下水流進行土壤絕對滲透性模擬。滲流模擬時將土樣相對的兩面作為出入邊界,假設流體在其余流動邊界及孔隙空間內(nèi)不產(chǎn)生滑移,且整個孔隙空間與外部隔離,流體不發(fā)生內(nèi)外交換。對出入邊界賦予壓力梯度,設定流體的流入及流出端。本文將邊界條件設置為入口壓力130 kPa,出口壓力100 kPa,水動力粘度為0.001 Pa·s,進行Z軸向的滲流模擬。在模擬時,對出入邊界上流體的流動速度進行積分可以得到流體流通體積,根據(jù)達西定律獲得試樣的絕對滲透率,計算式為

(1)

式中K——模擬得到的試樣絕對滲透率,μm2

μ——流體動力粘度,Pa·s

L——試樣在模擬流動方向的長度,cm

A——試樣垂直于流動方向的截面積,cm2

ΔP——出入邊界處的壓力差,kPa

Q——流體流通體積,cm3/s

為驗證單相滲流模擬結果的可靠性,需要將其與實際試驗結果進行對比,因此需要將絕對滲透率換算為絕對滲透系數(shù)(飽和導水率),換算關系為

(2)

式中ρ——流體密度,g/cm3

g——重力加速度,取10 m/s2

模擬結果采用樣本誤差比幾何平均數(shù)(GMER)、幾何標準偏差(GSDER)、均方根誤差(RMSE)進行評價[19]。

2 結果與分析

2.1 孔隙網(wǎng)絡模型參數(shù)

2.1.1孔隙半徑

土壤微觀結構是影響土壤滲透率的重要因素,其中孔徑分布可以詳細表征孔隙結構特征[20]。由圖3可知,不同土壤孔徑分布存在顯著差異。隨圍墾年限的增加,濱?？h土樣各土層的連通孔隙半徑均呈現(xiàn)增大的趨勢,較大連通孔隙對孔隙體積的貢獻率也相應提高。其中Y3處理除40～60 cm土層外,連通孔隙半徑多集中在100～500 μm之間。如東縣土樣中Y3處理相較其他處理,各土層中大尺寸連通孔隙含量增加。在Y1處理中各土層中對總孔隙體積貢獻率最大的孔隙集中分布在400～500 μm、300～400 μm、200～300 μm范圍。東臺市土樣在Y3處理中,0～40 cm土層土壤連通孔隙集中在100～400 μm,40～60 cm土層中連通孔隙多集中在100～300 μm,表明沿著深度方向,土壤連通孔隙半徑呈下降趨勢。

圖3 REV連通孔隙半徑分布及孔隙體積貢獻率

2.1.2配位數(shù)

孔隙空間中與單個獨立孔隙相連的喉道數(shù)量被稱為孔隙體配位數(shù)(CN)[21]。由圖4可知,各土樣孔隙配位數(shù)分布較為一致,均呈現(xiàn)較小配位數(shù)分布頻率較高,隨配位數(shù)增大其分布頻率顯著降低。濱海縣土樣Y1～Y3處理,各深度土壤孔隙配位數(shù)頻率分布曲線逐漸變緩,表明隨圍墾時間增加,土壤孔隙平均配位數(shù)有增加趨勢,孔隙局部連通性增強。同樣地,東臺市土樣在0～20 cm土層中,CN大于等于6時配位數(shù)頻率分布由大到小依次為Y3、Y2、Y1。不同土層孔隙配位數(shù)呈現(xiàn)的規(guī)律為當CN大于等于6時,0～20 cm土層的配位數(shù)頻率分布最大,占孔隙總數(shù)的10.91%;隨土層深度增加,配位數(shù)頻率分布降低,40～60 cm土層的配位數(shù)頻率分布最小,占孔隙總數(shù)的2.56%,孔隙空間的局部連通性隨土層深度增加而降低。

圖4 REV孔隙配位數(shù)分布

2.1.3喉道尺寸

喉道是相連孔隙之間的狹窄通道,喉道半徑表征流體遷移通道的大小,在一定程度上反映流體滲透的難易程度。隨圍墾時間的增加,較大半徑喉道占比有所提高(圖5)。同時在3種深度處理中,0～20 cm土層中半徑大于200 μm的喉道占比在全部深度處理中均為最高。東臺市0～40 cm土層相較于40～60 cm土層,土壤喉道半徑分布范圍更廣。如東縣土樣除Y1處理外,半徑大于200 μm的喉道占比隨土層深度增加而降低。

圖5 REV喉道半徑分布

圖6為各REV的喉道長度頻率分布及累積頻率分布。各土樣喉道長度頻率分布呈現(xiàn)單峰曲線。濱?？h土樣隨土層深度增加,喉道長度有增大的趨勢,各土層的平均喉道長度為964.08、991.81、1 206.63 μm。東臺市土樣隨圍墾時間的增加,各土層的平均喉道長度分別為1 485.29、1 635.97、1 000.09 μm。如東縣土樣在不同土層深度下,孔隙空間內(nèi)喉道長度大于1 000 μm的長喉道占比結果為:0～20 cm土層占比最小,占喉道總數(shù)的38.49%;40～60 cm土層占比最大,占喉道總數(shù)45.89%。

圖6 REV喉道長度分布

2.2 滲透性模擬

利用構建的孔隙網(wǎng)絡模型,模擬飽和單相水在REV中的滲流過程,得到各土樣的模擬滲透系數(shù)(飽和導水率)。由圖7a～7c可以發(fā)現(xiàn),絕對滲透系數(shù)(飽和導水率)的模擬值基本保留了實測值的數(shù)據(jù)規(guī)律特點,大致呈現(xiàn)出隨圍墾年限增加而升高的趨勢;而在土壤剖面上,則表現(xiàn)為隨土層深度增加而降低。同時對絕對滲透系數(shù)(飽和導水率)的模擬值進行評價分析(圖7d),結果顯示模擬值和實測值之間吻合度較高;GMER小于1,表明模擬結果低于實測值;具有較高絕對滲透系數(shù)(飽和導水率)的土樣,其模擬得到的結果誤差較大。

圖7 REV模擬絕對滲透系數(shù)(飽和導水率)和實測結果對比

2.3 土壤滲透系數(shù)的主要影響因素

選取部分孔隙網(wǎng)絡結構參數(shù)和絕對滲透系數(shù)模擬值進行灰色關聯(lián)分析,定量描述孔隙結構各參數(shù)對土壤滲透能力的影響。關聯(lián)度范圍為0～1,關聯(lián)度越大,因素間關系越密切[22-23]。計算選取的15個孔隙結構參數(shù)及其與絕對滲透系數(shù)的關聯(lián)度結果

如表1所示。各參數(shù)與絕對滲透系數(shù)的關聯(lián)度均在0.78以上。與土壤絕對滲透系數(shù)關聯(lián)度排序前5的參數(shù)由高到低分別為:最大喉道半徑、平均喉道半徑、孔隙平均配位數(shù)、孔隙率、平均孔隙半徑;喉道半徑參數(shù)對土壤絕對滲透系數(shù)的影響最為明顯。

表1 各參數(shù)與模擬絕對滲透系數(shù)(飽和導水率)的關聯(lián)度

3 討論

本文通過構建孔隙網(wǎng)絡模型,利用幾何體來表示孔隙空間,在一定程度上簡化了土壤結構,提高了模型的運行效率,減少了內(nèi)存存儲的計算壓力[24]。LI等[15]通過提取土壤孔隙網(wǎng)絡結構,分析了馬蘭黃土的孔隙度、孔徑分布及三維傾角等結構參數(shù),深入了解了黃土內(nèi)部微觀結構。WEI等[14]分析不同地層黃土孔隙網(wǎng)絡結構,發(fā)現(xiàn)淺層土壤具有網(wǎng)絡結構更加復雜、孔徑范圍更加寬泛等特征。在本研究中,通過網(wǎng)絡結構提取,可以發(fā)現(xiàn)不同土層深度和圍墾年限的土壤孔隙網(wǎng)絡差異顯著,淺層土壤孔隙擁有較大的連通孔隙尺寸范圍、較大平均配位數(shù)、較大的喉道半徑以及較小的喉道長度;同樣地,圍墾時間較長土壤的連通孔隙分布更為廣泛、平均配位數(shù)及喉道半徑均更大,喉道長度也相對較小,表明農(nóng)業(yè)圍墾活動在一定程度上改善了沿海圍墾區(qū)土壤結構。對比WEI等[14]的研究,發(fā)現(xiàn)本試驗中土樣孔隙網(wǎng)絡的平均孔喉半徑均大于不同地層黃土,這可能歸因于CT掃描的分辨率和選取REV尺寸的差異;但各土樣的平均配位數(shù)卻均遠小于黃土,這是因為相較于黃土,濱海土壤含有較高的交換性鈉離子,會對土壤膠體顆粒產(chǎn)生分散作用,導致土壤結構較為破碎、土壤連通性較差。

利用提取的孔隙網(wǎng)絡進行單相水滲流模擬,各土樣模擬絕對滲透系數(shù)和試驗測得的飽和導水率之間呈現(xiàn)出相對一致的變化規(guī)律,模擬值略低于試驗測定結果。原因可能為:①孔隙結構提取過程中受到圖像分辨率的限制,小于分辨率的連通孔隙不能被識別,因此也不能參與模擬過程。②在試驗過程中,土樣中部分孔隙會在水流作用下相互連通,形成新的孔隙通道,進而在模擬和試驗結果中產(chǎn)生差異[14-15]。

孔隙結構對多孔介質(zhì)中流體體積、分布及運移有重要影響,是影響流體滲流的直接因素。孔隙網(wǎng)絡參數(shù)能夠幫助解釋水分在土壤等多孔介質(zhì)的傳輸特征。過去的研究也表明,在孔隙尺度上,孔隙尺寸分布及孔隙拓撲關系是影響土壤水力特性和溶質(zhì)運移的兩個主控因素[25]。灰色關聯(lián)度分析顯示,最大喉道半徑與模擬絕對滲透系數(shù)的關聯(lián)度最高,其次是平均喉道半徑。同時本研究中平均喉道半徑、最大喉道半徑、最小喉道半徑的關聯(lián)度之和,大于平均孔隙半徑、最大孔隙半徑、最小孔隙半徑的關聯(lián)度之和。這一結果表明孔隙喉道半徑對土壤滲透性的影響作用最大。這也驗證了WEI等[14]的研究結論,其認為土壤中喉道半徑對控制水分流動起到主要作用,同時大喉道的影響作用高于多個小喉道的疊加效果。閆國亮等[12]發(fā)現(xiàn)喉道半徑對儲層巖石滲透率的影響大于孔隙半徑的影響。SONG等[24]也認為孔隙喉道結構是影響油儲層滲透能力的主要決定因素,較大的喉道以及良好的連通性可以為儲層內(nèi)致密油移動提供通道。配位數(shù)則是表征孔隙網(wǎng)絡連通性的重要參數(shù),而連通性會直接影響到土壤的滲透能力。因此通常認為平均配位數(shù)越大,土壤滲透系數(shù)越大。與此同時,水分傳輸移動的難易程度一般被認為與單位土壤中孔隙的體積密切相關,表現(xiàn)為孔隙率越大,土壤的水分運輸能力越強[14]。

4 結論

(1)圍墾年限和土層深度均會對土壤孔隙分布及其拓撲結構產(chǎn)生影響。隨圍墾年限的增加,土壤連通孔隙半徑范圍、平均配位數(shù)、喉道半徑增加,喉道長度減小;淺層土壤連通孔隙分布較為廣泛、平均配位數(shù)及喉道半徑均較大,喉道長度也相對較小。

(2)基于孔隙網(wǎng)絡模擬得到的絕對滲透系數(shù)和試驗測定的飽和導水率結果吻合度高,呈現(xiàn)一致的變化趨勢,孔隙網(wǎng)絡模型在表征土壤孔隙結構方面具有較好適用性和可靠性。

(3)最大喉道半徑與模擬絕對滲透系數(shù)的關聯(lián)度最高,表明喉道半徑對土壤滲透性能的影響作用最大。