包氣帶水分遷移離心模擬的可行性探討

(東華理工大學 核資源與環(huán)境國家重點實驗室，南昌 330013)

1 研究背景

包氣帶作為聯(lián)系“四水”(大氣水、地表水、包氣帶水和地下水)動態(tài)水循環(huán)系統(tǒng)的紐帶，是指從地下潛水面向上延伸至地面的地帶[1-2]，是大氣水和地表水同地下水發(fā)生聯(lián)系并進行水分交換的地帶，由于污染物大多是隨水擴散發(fā)生遷移，因此它也是地表污染物向下遷移污染地下水的通道[3]。近年來，人類活動的影響及各種污染物的排放，加劇了土壤的污染，土壤退化、地下水污染、荒漠化、地面沉降等一系列環(huán)境問題都與包氣帶土壤水分有密切聯(lián)系[4]。由于包氣帶在水文循環(huán)中處于重要位置，與降水入滲、地表蒸發(fā)、地下水補給、潛水蒸發(fā)等水文過程密切相關，同時，包氣帶水分運動是污染物運移及水、汽、熱運動的主要驅(qū)動力，因此，對包氣帶中水分的遷移進行研究，探明水分運動規(guī)律，對于研究“四水”轉(zhuǎn)換、地下水來源、水資源開發(fā)利用、生態(tài)保護及人類的生產(chǎn)生活都具有十分重要的意義[2,4-5]。

離心模擬是將按比例縮小尺寸的原型放在離心機內(nèi)，通過離心加速還原原型中應力水平并開展模型模擬的一種技術(shù)[6-7]。離心模擬通過加速水分和溶質(zhì)在土壤中的遷移速度，能夠大大縮短實驗時間[8]。相比于常規(guī)物理模擬方法，離心模擬具有快速搜集實驗數(shù)據(jù)、應力水平保持和原型一致、預測全尺度的遷移行為和可以開展二、三維模擬等優(yōu)點[6,9]。離心模擬技術(shù)已經(jīng)廣泛應用于包氣帶溶質(zhì)遷移的研究[8,10-13]，在低滲透土壤長期污染物遷移及獲取數(shù)據(jù)驗證模型等方面具有重要意義。然而，對于離心模擬技術(shù)能否應用于包氣帶水分遷移尚存爭議[6,14]。Goforth[15]和Poulose等[16]對該技術(shù)應用于包氣帶水分遷移表示懷疑，他們認為，在飽和度較低的土中，離心機不能正確模擬基質(zhì)勢作用下的流體流動。Mitchell[17-18]認為離心模擬能夠正確描繪原型污染物遷移現(xiàn)象，包括非飽和細砂土中基質(zhì)吸力引起的遷移。

水力傳導系數(shù)與含水率的關系K(θ)及土水特征曲線SWRC(Soil Water Retention Curve)對于理解包氣帶水分的存儲及遷移具有重要意義，常規(guī)方法確定K(θ)和SWRC存在許多挑戰(zhàn)，最主要的挑戰(zhàn)是實驗耗時過長，因此大多數(shù)包氣帶水分遷移研究都是依靠經(jīng)驗公式或理論模型?，F(xiàn)有研究表明，純粹依靠經(jīng)驗關系或理論模型并不可靠，通過實驗獲取非飽和土水力特性十分重要。Singh等[19]利用離心機測定了多種不飽和土的水力傳導系數(shù)及其特征曲線，發(fā)現(xiàn)與理論值很接近。

目前，離心模擬技術(shù)在包氣帶溶質(zhì)遷移的應用中都默認了該技術(shù)能夠成功模擬包氣帶水分遷移。上述爭議無疑會阻礙離心模擬技術(shù)的推廣和發(fā)展，然而目前并沒有相應的研究來解決這一爭議。本文通過建立離心場內(nèi)包氣帶水分遷移的一維理論模型和數(shù)值模型，選用4套文獻數(shù)據(jù)對數(shù)值模型進行驗證，在此基礎上探討離心模擬技術(shù)應用于包氣帶水分遷移研究的可行性，期望從理論上對該爭議的澄清進行探討。

2 包氣帶水分遷移離心模型

離心模型中非飽和水流是受水勢梯度驅(qū)動的，其控制方程為Richards方程。由于離心模型處于離心力場內(nèi)，故離心條件下包氣帶水分遷移的控制方程形式不同于自然重力下常規(guī)的情況。

2.1 理論模型

2.1.1 水勢與達西定律

圖1是一維離心模型中非飽和水流的原理分析圖，其中rt表示模型頂部半徑；rb表示模型底部半徑；z表示當前位置離參照面的距離；L表示離心模型的高度。驅(qū)動水流運動的水勢主要有水的位置勢能、動能和基質(zhì)勢能3部分。

圖1 水分遷移一維離心模型分析圖Fig.1 Analysis of one-dimensional water migration centrifugal model

離心模型中水分的滲流速度一般來說很小，水的動能可忽略不計，離心模型中單位質(zhì)量水的水勢能為

Φ=Pe+ψm/ρ。

(1)

式中：Φ為水的勢能；Pe為離心狀態(tài)的水的位置勢能；ψ是基質(zhì)吸力；ρ是水的密度；下標m代表離心模型。

離心狀態(tài)下和常規(guī)情況下水的勢能的主要差別在于位置勢能的計算不同。離心加速度沿徑向存在一個分布趨勢而不是某個固定值，其大小是半徑和角速度的函數(shù)(a=rω2)。離心狀態(tài)下水的位置勢能等于將單位質(zhì)量的水從參照面移至當前位置需克服離心力所做的功，將模型底部設為參照面，則半徑r處單位質(zhì)量水的離心位置勢能為

式中：a為離心加速度；ω為角速度。

聯(lián)合式(1)和式(2)可得一維離心模型非飽和水的勢能為

(3)

達西定律描述了水流運動和水勢梯度之間的關系，離心情況下達西定律的適用性已被許多研究證實[8,13]，其表達式為

(4)

(5)

圖2 一維離心模擬的節(jié)點剖分 Fig.2 Node dissection of the one-dimensional centrifugal modeling

2.1.2 Richards方程

考慮離心模型中某一個控制體積單元(圖1)，根據(jù)連續(xù)性原理有(t為離心模擬時間；θm是離心狀態(tài)下的含水率)

(6)

將式(5)代入式(6)可得

式(7)是包氣帶水流的控制方程，它刻畫了離心場內(nèi)水流流經(jīng)非飽和土壤的過程規(guī)律。

2.2 數(shù)值模型

2.2.1 模型數(shù)值化

用n個節(jié)點將一維離心模型分成n個有限差分網(wǎng)格(圖2)，圖中第n+1個節(jié)點是虛擬節(jié)點，是為了在求解中刻畫離心模型頂部邊界條件所設。對于模型中第i個節(jié)點，利用改進的Picard循環(huán)[20]線性化式(7)中的非線性項，則等號后面項的有限差分近似為

(8)

對時間項歐拉形式向后差分，并嵌入Picard循環(huán)(Δt為時間步長)得

(9)

定義土壤的特殊水容量為

(11)

聯(lián)合式(9)—式(11)，則式(7)有限差分為

聯(lián)立式(8)和式(12)，將未知項放在左邊，已知項放在右邊，則第j+1個時間步、τ+1個Picard循環(huán)中基質(zhì)吸力可通過式(13)求得,即

(13)

將式(13)應用于離心模型上所有節(jié)點，得到一個方程組，可用追趕法求解。對于邊界節(jié)點，式(13)需要調(diào)整以反映不同的邊界條件。

2.2.2 邊界條件和物料守恒

根據(jù)文獻[20]，數(shù)值模擬時底部的自由流邊界可處理成定水頭邊界。因此，本文考慮的邊界條件為：模型頂部提供一個穩(wěn)定通量的水流，底部是一個定基質(zhì)吸力邊界，即

(14)

式中：qb是頂部的水流通量；ψb是底部穩(wěn)定的基質(zhì)吸力。

為了了解一維離心模型中非飽和水流的物料守恒情況，需要知道上下邊界的水流通量和模型中含水率的變化。模型底部的出流水通量可以利用泰勒級數(shù)推導而出，即

(16)

表1 例1—例4數(shù)值模擬參數(shù)取值Table 1 Parameter values of examples 1-4 in numerical modeling

注：?表示文獻中沒有直接給出而是從文獻的圖表中獲取的參數(shù)；表示原文獻該實驗重復了2次，表中所給數(shù)據(jù)為2次實驗的平均值

2.2.3 土壤水力學特性

求解Richards方程要知道水力傳導系數(shù)和基質(zhì)吸力、含水率和基質(zhì)吸力之間的關系[21]，即滲透率函數(shù)K(θ)和土水特征曲線SWRC。常用Van Genuchten模型[22]描述SWRC，即

(18)

式中：Θ是有效飽和度；α，nv，mv=1-(1/nv)是Van Genuchten模型參數(shù)；θs和θr分別為飽和含水率和殘余含水率；h是壓力水頭，可用h=ψ/(ρg)實現(xiàn)和基質(zhì)吸力的轉(zhuǎn)換。

K(θ)有許多種形式，其中Gardner模型[23]和Mualem模型[24]是最常用的，即：

(19)

(20)

式中：Ks是土壤的飽和水力傳導系數(shù)；ψo和β是Gardner模型參數(shù)；l是孔隙連通性參數(shù)，對大部分土壤來說為0.5[24]。

3 包氣帶水分遷移離心模擬可行性探討

3.1 文獻數(shù)據(jù)的驗證

非飽和水分遷移離心數(shù)值模型用4套已發(fā)表的文獻數(shù)據(jù)(表1)進行驗證，參數(shù)取值均從文獻中獲取(有些通過圖表獲取實驗數(shù)據(jù)，再利用Solver求解器間接擬合得到)。

3.1.1 例1：ω=0時基質(zhì)吸力主導的非飽和水流

當沒有角速度時，離心加速度為0，非飽和水的位置勢能處處相等，水流運動受基質(zhì)吸力梯度主導。本套數(shù)據(jù)來自Kirkham[25]，使用一根總長25 cm由透明有機玻璃制作的水平土柱開展實驗，整根土柱填滿紅色Ferrosol土，填充重度為0.155 g/cm3。土柱一端連接一個馬氏瓶，另外一端保持開放暴露在空氣中。由于實驗是破壞性取樣，為了獲得一個時間系列的數(shù)據(jù)共做了6次實驗，初始和邊界條件不能完全相同，本文采用6次實驗初始和邊界條件的平均值作為參數(shù)輸入，并將連接馬氏瓶的一端處理成固定基質(zhì)吸力為0 kPa的邊界條件，另外一端處理成自由流邊界。模擬結(jié)果和實驗數(shù)據(jù)對比(圖3)說明兩者匹配較好。

圖3 水平土柱實驗基質(zhì)吸力梯度主導的非飽和水流實驗數(shù)據(jù)[25]和模擬結(jié)果的對比Fig.3 Comparison between the experimental data ofunsaturated flow dominated by matrix suction gradient andthe simulated results in horizontal soil column experiment

3.1.2 例2：穩(wěn)態(tài)流離心模型剖面水分曲線

本套數(shù)據(jù)來自Nimmo等[26]，用來驗證本文提出的理論模型在預測穩(wěn)態(tài)模型剖面上水分分布曲線方面的能力。實驗過程中土柱部分內(nèi)填充Oakley砂，預先用脫氣的0.01N CaSO4和0.01N CaSeO4混合液飽和，在100 rad/s(對應離心加速度約194g)的角速度下離心。離心過程中頂部提供穩(wěn)定的水流通量，底部維持固定基質(zhì)吸力為-10 kPa的邊界條件。離心穩(wěn)定至沒有水分出流后，土樣被切成1.25 mm厚的片測定不同部位的含水率。此處利用數(shù)值模型根據(jù)其實驗狀態(tài)來預測穩(wěn)態(tài)剖面的水分分布情況，模擬結(jié)果和實驗數(shù)據(jù)對比如圖4(a)，可以看到兩者相當吻合。圖4(b)是預測的基質(zhì)吸力分布曲線，從中可以看出當穩(wěn)態(tài)建立后，靠近離心模型頂部的水分及基質(zhì)吸力分布趨近均一，重力水平N值越大現(xiàn)象越明顯[20]。

圖4 穩(wěn)態(tài)剖面曲線Fig.4 Numerical simulation of the steady-state profile curve

3.1.3 例3：暫態(tài)非飽和水流

本套數(shù)據(jù)來源于Nakajima等[27]，他們使用一臺半徑2 m、安裝有基質(zhì)勢傳感器且能夠測量累計出流量的土工離心機開展實驗，擬利用一步出流法估算非飽和土的特性參數(shù)。離心過程中土樣頂部水流通量為0 m/s，底部保持和水面的接觸(可視作固定基質(zhì)吸力為0 kPa的邊界條件)。實驗使用了10，20，40g的離心條件，但調(diào)整了模型的尺寸，保證了模擬的是同一個原型。Nakajima等[27]利用獲得的實驗數(shù)據(jù)估算了非飽和土特性參數(shù)，本文利用他們估算的實驗參數(shù)逆向模擬暫態(tài)的非飽和水分遷移過程。圖5對比了模擬結(jié)果和實驗數(shù)據(jù)，兩者有著較為滿意的吻合度。40g條件下匹配較差的原因有2個：不理想的邊界條件控制和基質(zhì)吸力傳感器(張力計)性能限制。從圖5可以看到在原型時間為80 min時，土樣底部的基質(zhì)吸力水頭為正值，表明土樣底部排水不暢甚至發(fā)生了積水現(xiàn)象。此外，重力水平N較大的情況下，會引起快速的孔隙水壓變化，會對張力計的陶瓷頭產(chǎn)生破壞從而引起讀數(shù)不準。

圖5 數(shù)值模擬結(jié)果和實驗數(shù)據(jù)對比Fig.5 Comparison between numerical simulation resultsand experimental data

3.1.4 例4：多重轉(zhuǎn)速離心實驗非飽和水流

表2 Run 1的離心轉(zhuǎn)速改變時間序列Table 2 Time series of the centrifugal speed changes in Run 1

土樣底部的固定基質(zhì)吸力在一次多重轉(zhuǎn)速離心實驗中并不是全程不變的，每改變一次離心轉(zhuǎn)速，該值就相應改變一次。本文假設多孔陶瓷片在離心過程中每一處的水勢相等，土樣底部的固定基質(zhì)吸力可通過式(21)計算。

(21)

式中H為陶瓷片的厚度。利用數(shù)值模型預測了Run 1實驗過程中的水分遷移變化并和實驗數(shù)據(jù)比較(圖6(a))，模擬結(jié)果和實驗數(shù)據(jù)很吻合。此外，本文還預測了在每個轉(zhuǎn)速改變時間的節(jié)點上，土樣的基質(zhì)吸力剖面分布曲線(圖6(b))，和等[28]給出的完全一致，進一步說明了理論模型的正確性。

以上4個不同場景的實驗可以充分驗證理論及數(shù)值模型的正確性，但仍不足表明數(shù)值模型能夠正確預測離心非飽和流。以例1為例，結(jié)合式(17)考察數(shù)值模型的物料守恒特性，結(jié)果如圖7所示。本文將物料守恒誤差在5%內(nèi)視為優(yōu)秀，根據(jù)圖7可知當時間步長和空間步長均足夠小時，數(shù)值模型的物料守恒特性優(yōu)秀。此外，數(shù)值模型對空間步長比時間步長更敏感：當時間步長從1 s跨越到60 s時，96 min時刻的非飽和水分遷移的物料誤差均在5%以內(nèi)；當空間步分別長為2.5 mm和5 mm時，物料誤差分別達到并超過了10%和60%。因此，當利用數(shù)值模型時，建議使用適當?shù)臅r間步長提高計算效率，使用較小的空間步長保證模擬的質(zhì)量。要強調(diào)的是情景1不是處在離心場內(nèi)，而在離心場下應當使用更小的空間步長，尤其是靠近底部出流邊界處的空間步長更應如此[28]。

圖6 Run 1實驗模型預測結(jié)果和實驗數(shù)據(jù)對比及土樣基質(zhì)吸力剖面分布曲線Fig.6 Comparison between experimental model prediction results and experimental data in Run 1

圖7 以例1為例考察數(shù)值模型的物料守恒特性Fig.7 Test of materials conservation in the numerical model (example 1)

3.2 可行性爭議與探討

如前所述，Goforth[15]基于①非飽和狀態(tài)下滲流速度的離心相似比可能不成立；②常規(guī)非飽和水流中，基質(zhì)吸力梯度可能是重力勢能梯度的10～1 000倍，認為將離心機應用于包氣帶水分遷移不太合適。Poulose等[16]通過實驗觀察認同Goforth[15]的結(jié)論，認為離心模擬只適合飽和水分及溶質(zhì)遷移的研究。針對這一爭議，本文按照以下步驟進行探討：首先，認為“非飽和流的滲流速度、尺寸和時間的相似比分別為N，1/N，1/N2”是正確的；其次，為了判斷離心模擬能否正確重現(xiàn)原型現(xiàn)象，利用Hydrus 1-D軟件包直接對全尺度原型中非飽和水分遷移過程進行模擬預測；最后，將兩者的模擬結(jié)果進行比較從而進一步分析判斷。

此處的原型為一根均質(zhì)裝填低滲透性非飽和土土柱的水分入滲過程。土柱的初始含水率設為0.22 cm3/cm3，入滲水從頂部注入并保持通量恒定，土柱底部保持和水面接觸，具體的輸入?yún)?shù)見表3。

數(shù)值模型的模擬結(jié)果和Hydrus軟件包直接模擬原型的結(jié)果對比如圖8所示，可以看出兩者幾乎完全一致，表明離心模擬非飽和水流的相似比是正確的。細心觀察圖8，可看到在靠近底部時，模型和原型的數(shù)據(jù)匹配存在一定的偏差，離心模型的水分遷移要稍滯后于原型，這是由離心加速度的不均勻分布所導致的[3]。模型2/3高度處的相似比為N，而柱頂至2/3高度處離心加速度略小于Ng(相似比N)。另外，用同一離心設備模擬同一原型時，使用更高的離心加速度能夠減輕這種現(xiàn)象。因為更大的N需要更小尺寸的離心模型，離心加速度分布更加均勻。

圖8離心機用于非飽和水流可行性的數(shù)值檢驗
Fig.8Feasibilityofapplyingcentrifugeinunsaturatedwaterflowbynumericaltests

針對Goforth[15]給出的2個原因，本文從以下2方面進行探討。

表3 數(shù)值模擬的輸入?yún)?shù)Table 3 Input parameters of the numerical modeling

一方面，Goforth[15]忽略了能夠滿足滲流速率正比于重力水平的另一種情況，即基質(zhì)吸力梯度(?ψm/?r)正比于重力水平。從圖8可知，離心模型和原型對應位置處的基質(zhì)吸力大小相等，離心模型的尺寸和原型相比縮小了1/N，所以離心模型中2點之間的基質(zhì)吸力梯度大小應該是原型中對應2點間梯度的N倍，即基質(zhì)吸力梯度和N成正比。故滲流速度的相似比為N是正確的，離心模擬技術(shù)可以正確反映原型中的非飽和水分遷移。

另一方面，離心模型和原型中位置勢能梯度占總水勢梯度的比分別為：

(22)

(23)

由于基質(zhì)吸力梯度正比于重力水平N，所以ηm和ηp相等，這意味著離心不會強化位置勢能梯度在非飽和帶水分遷移過程中的重要性。

鑒于以上原因，本文認為Goforth[15]思考不嚴謹，其結(jié)論是不正確的。本文中的理論和數(shù)值模型都假設土的特性參數(shù)不隨離心加速度的改變而改變，基于這一前提，可以說理論上離心模擬能用于包氣帶水分遷移的研究。但是在實際中存在許多因素能導致不確定性，如不理想的邊界條件控制、壓實現(xiàn)象的出現(xiàn)、安裝的傳感器性能、離心模型制作、加速度分布不均等，在進行包氣帶水分遷移離心模擬研究時，應該重視解決這些問題確保得到準確的結(jié)果。

3.3 理論分析與數(shù)值模擬的不足

本文針對離心模擬技術(shù)應用于包氣帶水分遷移研究的可行性而進行的理論分析和數(shù)值模擬依然存在其不足之處，主要表現(xiàn)為：①理論模型中忽略了水的動能，這對于一維離心模型非飽和水的勢能的推導會產(chǎn)生影響，進而會對后續(xù)的數(shù)值模擬的結(jié)果產(chǎn)生影響，這種影響使得理論分析并沒有完全描述包氣帶水分遷移的過程；②數(shù)值模擬過程中求解Richards方程要知道水力傳導系數(shù)和基質(zhì)吸力、含水率和基質(zhì)吸力之間的關系，即滲透率函數(shù)K(θ)和土水特征曲線SWRC，而本文用于描述K(θ)和SWRC的模型是常用的Gardner模型和Mualem模型及Van Genuchten模型，此處并沒有論證其它描述K(θ)和SWRC的模型是否可用于本文，這一點值得進行相關的探討；③數(shù)值模擬部分驗證了從文獻中獲取的4套數(shù)據(jù)，雖然這4套數(shù)據(jù)對應的4種情景足夠?qū)Ρ疚牡淖h題進行探討和驗證，但在后續(xù)的進一步研究中，數(shù)值模擬部分可以將驗證的情景進行擴展，以便更加深入地進行相關研究的探討，本文作者在后續(xù)的研究中會對此進行考慮。

4 結(jié) 論

本文建立了離心場內(nèi)一維包氣帶水分遷移的理論和數(shù)值模型，選用4套文獻數(shù)據(jù)對數(shù)值模型進行了驗證，在此基礎上探討了將離心機用于包氣帶水分遷移研究的可行性。研究結(jié)論如下：

(1)數(shù)值模型能夠很好地重現(xiàn)文獻中包氣帶水分遷移過程，與實驗數(shù)據(jù)匹配較好，說明了理論模型的正確性和數(shù)值模型的準確可靠性。

(2)分析物料守恒誤差發(fā)現(xiàn)，當時間步長和空間步長取值都足夠小時，數(shù)值模擬結(jié)果的質(zhì)量能夠保證為優(yōu)秀，在實際的數(shù)值模擬中建議采用小的空間步長(確保模擬結(jié)果的質(zhì)量)和較大的時間步長(提高計算效率)。

(3)離心模擬技術(shù)應用于包氣帶水分遷移是可行的。Goforth[15]的結(jié)論不嚴謹，因為基質(zhì)吸力梯度和重力水平N成正比，離心不會強化位置勢能梯度在驅(qū)動非飽和水分遷移過程中的重要性，滲流速度、尺寸和時間的相似比分別為N，1/N，1/N2是正確的。

(4)離心模擬時加速度分布不均會導致水分模型底部區(qū)域的遷移滯后于原型，在用同一設備模擬同一原型時，使用更高的離心加速度能夠減輕這種現(xiàn)象。

(5)離心模擬時一般假設土的特性參數(shù)不隨離心加速度的改變而改變，然而，實際中存在許多因素導致實驗數(shù)據(jù)和理論產(chǎn)生較大偏差，包括不理想的邊界條件控制、壓實現(xiàn)象的出現(xiàn)、安裝的傳感器性能、離心模型制備、加速度分布不均等。