孫露露，楊衛(wèi)波

(揚州大學水利與能源動力工程學院，江蘇揚州 225127)

溫度梯度作用下非飽和土壤中溶質遷移過程的數(shù)值模擬

孫露露，楊衛(wèi)波*

(揚州大學水利與能源動力工程學院，江蘇揚州 225127)

根據非飽和土壤非等溫Richard方程、傳熱與溶質傳遞方程，建立溶質遷移模型，研究溫度梯度作用下入滲速度對土壤溫度變化及溶質遷移規(guī)律的影響.數(shù)值模擬的結果表明:溶質遷移過程中，土壤溫度與溶質的物質的量濃度均隨時間延續(xù)而逐漸升高，隨深度的加深而逐漸降低;入滲速度的增加可促進土壤中熱量的擴散和對流，增強溶質的對流和彌散作用.

非飽和土壤;溫度梯度;滲流;溶質遷移;數(shù)值模擬

近年來，溶質在土壤中的遷移機制已成為當前土壤環(huán)境控制領域的一個研究熱點.非飽和土壤是由固、液、氣三相組成的多孔介質，溶質在沿非飽和土壤的水流運動過程中會導致土壤污染.Viotti等［1］建立了一維非飽和土壤中溶質遷移方程，研究了不同輸入參數(shù)對溶質濃度分布的影響;Kumar等［2］通過無網格Galerkin(element free galerkin，EFG)法分析了二維非飽和多孔介質中溶質遷移問題，并與有限元法進行了算例比較;Javadi等［3］對非飽和土壤中水、氣和溶質的遷移進行了案例分析;Ghasemzadeh［4］研究了變形多孔介質中熱量、水流和溶質的遷移過程;Winiarski等［5］實驗研究了非飽和土壤的非均質性對多相流和溶質遷移特性的影響;眭素剛等［6］開展了非飽和土壤彌散、滲流特性試驗，并模擬分析了分層土壤中溶質的濃度分布;左自波等［7］研究了不同降雨條件下非飽和帶孔隙水壓力、含水量變化及溶質的遷移規(guī)律.鑒于溶質類型的多樣性，本文以某一特定濃度的無機液體A為代表，在忽略吸附、衰減及化學性質的條件下，根據多孔介質傳熱、傳質及滲流理論，建立非飽和土壤中溶質遷移的數(shù)學模型，研究入滲速度對非飽和土壤溫度及溶質的物質的量濃度c(A)分布規(guī)律的影響，以及滲流作用下非飽和土壤的傳熱傳質規(guī)律.

1 計算模型

1.1 物理模型

為模擬滲流影響下土壤溫度和溶質濃度的時空變化規(guī)律，以非飽和土壤為研究對象，建立溶質遷移模型，進行滲流場、溫度場和溶質傳遞場三場耦合的二維數(shù)值模擬，研究溫度及溶質的物質的量濃度的分布規(guī)律，模擬對象結構如圖1所示.圖1中ABCD區(qū)域為非飽和土壤，模型尺寸為AB=1 m，AD=0.5 m;E(0.5，0.45)為溶質源.溶質遷移在豎直方向上分別取5個觀測點，1#(0.5，0.425)，2#(0.5，0.4)，3#(0.5，0.35)，4#(0.5，0.25)，5#(0.5，0.1).

圖1 非飽和土壤中溶質遷移模型示意圖Fig.1 Schematic diagram of physicalmodelw ith solute transport in unsaturated soils

1.2 數(shù)學模型

由于非飽和土壤中滲流、傳熱及溶質傳遞是一個多因素互相耦合的復雜過程，為便于模型建立與求解，現(xiàn)作如下簡化假設:①土壤水不可壓縮;②土壤土質均勻，各向同性;③不存在氣體、固體和其他補給及化學反應來源;④水體的密度、導熱性和比熱容恒定;⑤忽略溫度變化對土壤熱物理性質的影響.

1)滲流控制方程.非飽和土壤滲流滿足Darcy定律［8］:

式中v為滲流速度;k為滲透系數(shù);H為水力水頭.

基于質量守恒原理，非飽和土壤二維滲流的非等溫Richard方程［9］為

式中θ為土壤體積含水率;τ為時間;kx，ky分別為土壤x與y方向的等溫滲透系數(shù);kTx，kTy分別為土壤x與y方向的非等溫滲透系數(shù);T為溫度.

非飽和土壤中體積含水率和土壤水吸力呈非線性關系，可由土壤物理學領域最為普遍使用的van Genuchten土—水特征曲線通用方程［10］

表示，式中θ為土壤體積含水率;θr為殘余體積含水率;θs為飽和體積含水率;ψ為土壤水吸力;n，m，α為模型參數(shù)，其中n為土—水特征曲線指數(shù)，與曲線的陡緩程度有關且滿足n=(1－m)－1.

將van Genuchten土—水特征曲線與土壤滲透系數(shù)Mualem模型［11］聯(lián)合，可得土壤滲透系數(shù)的解析表達式:

式中Se=(θ－θr)/(θs－θr)為有效飽和度;kr=k/ks為相對滲透系數(shù)，其中ks為飽和滲透系數(shù);m為模型參數(shù).

2)傳熱控制方程.當土壤中的滲流為穩(wěn)定流時，熱傳輸方程［8］330為

式中cs為土壤比熱容;λ為土壤熱導率;cw為土壤中流體(水)的比熱容;u為滲透速度.

3)溶質傳遞方程.基于溶質的質量守恒定律，可得非飽和土壤的二維溶質傳遞方程［8］231:

4)定解條件的確定.

①滲流問題.對于非飽和土壤模型，其頂部CD為速度入口邊界條件;AB為常水頭邊界條件，水頭為0;其余邊界為無流動邊界.非飽和土壤ABCD區(qū)域，水頭壓力初始值為負數(shù).

②傳熱問題.非飽和土壤模型的AB和CD為溫度邊界條件，分別為300 K和308 K，其余邊界為周期性熱條件邊界.非飽和土壤ABCD區(qū)域初始溫度為邊界條件下穩(wěn)定的溫度場，如圖2所示.

③溶質傳遞問題.非飽和土壤模型中點E為溶質源，通量邊界條件為 0.5 mol·(m2·s)－1，其余為無通量邊界條件.整個土壤模型區(qū)域水中溶質的初始濃度為0.

圖2 非飽和土壤溶質遷移前穩(wěn)態(tài)溫度場(K)Fig.2 Tem perature distribution before solute transport in unsaturated soils

1.3 計算參數(shù)及數(shù)值計算

采用Comsol Multiphysics多物理場耦合分析軟件中多孔介質地下水流模塊、傳熱模塊及物質傳遞模塊，進行3個物理場的耦合模擬計算.非飽和土壤滲流計算參數(shù)由土壤水特性軟件Spaw得到:滲透系數(shù)k為3.5 ×10－6m·s－1，飽和含水率 θs為0.43，殘余含水率 θr為 0.08，土壤密度 ρs為 1 400 kg·m－3，關系常數(shù)α，n分別為1.74和1.38;傳熱計算參數(shù):土壤熱導率 λs為 2 W·(m·K)－1，土壤比熱容 cs為 2 000 J·(kg·K)－1，水的密度 ρw為1 000 kg·m－3，水的熱導率 λw為0.54W·(m·K)－1，水的比熱容 cw為4 200 J·(kg·K)－1;溶質傳遞計算參數(shù):分子擴散系數(shù)D為8.25×10－7m2·s－1，水平彌散度αr為0.001 m，垂直彌散度αz為0.002 m.

本文采用三角形網格對非飽和土壤進行自由網格劃分，近溶質源區(qū)域網格較密，并進行了網格精度檢驗.由于基本假設中忽略了溫度變化對土壤熱物理性質的影響，且水體的物性參數(shù)恒定，所以實際大氣溫度波動對溶質遷移的影響較小，可忽略不計.恒定的環(huán)境溫度條件下，土壤具有穩(wěn)定且呈線性分布的溫度場，形成滲流條件下溶質遷移模型的初始溫度梯度，見圖2.

2 計算結果及分析

2.1 不同深度處土壤溫度及溶質的物質的量濃度分布

對入滲速度為2×10－6m·s－1時非飽和土壤中溶質遷移模型進行數(shù)值計算，分別得到第25 h時刻非飽和土壤溫度及溶質的物質的量濃度分布云圖、不同深度處非飽和土壤溫度及溶質的物質的量濃度變化曲線，見圖3～6.

圖3 第25 h時刻土壤溫度(K)分布云圖Fig.3 Tem perature distribution cloud graph at 25 h

圖4 第25 h時刻溶質的物質的量濃度分布云圖Fig.4 Concentration distribution cloud graph at 25 h

由圖3可觀察得到，土壤溫度由上至下逐步降低，與土壤初始時的溫度相比有了較大的提升.這主要是由于入滲雨水溫度高于土壤本身的初始溫度，隨著雨水滲流進入土壤，與土壤固體骨架及水分進行導熱和對流換熱，從而導致溫度逐漸升高.由圖4可以觀察得到，靠近溶質源區(qū)域溶質的c(A)越高，遠離溶質源區(qū)域c(A)越低，且溶質沿豎直方向的遷移較水平方向更大.這是由于溶質在土壤中的遷移與溶質彌散度、水流在土壤中的滲透速度相關，且豎直方向的溶質彌散度及水流的滲透速度均大于水平方向.

圖5 不同深度處土壤溫度變化曲線Fig.5 Variations of soil tem perature at different depths

圖6 不同深度處溶質的物質的量濃度變化曲線Fig.6 Variations of concentration at different depths

由圖5可觀察得到，1#處土壤溫度最高，5#處最低，且整體隨著時間增加而升高.在0～5 h期間，土壤溫度有一定波動，越深處波動越明顯.這是由于滲流前期土壤逐漸飽和，下層土壤含水率增大導致土壤熱導率減小，這與上層土壤熱導率差異較大;滲流后期土壤基本飽和，豎直方向熱導率差異逐漸縮小.由圖6可觀察得到，溶質的c(A)隨時間變化而不斷增大，越接近溶質源區(qū)域溶質的c(A)越高且增幅更明顯.這是由于越接近土壤地表，溶質的c(A)變化受滲流影響越顯著，且隨著深度增加影響越小.

2.2 不同入滲速度時土壤溫度及溶質的物質的量濃度分布

為分析入滲速度對非飽和土壤溫度及溶質的物質的量濃度分布特性的影響，對入滲速度分別為1×10－6，2×10－6，3×10－6m·s－1時溶質遷移模型進行數(shù)值計算，分別得到不同入滲速度時非飽和土壤中1#處的溫度及溶質的物質的量濃度變化曲線，見圖7～8.

圖7 不同入滲速度時1#處土壤溫度變化曲線Fig.7 Variations of soil tem perature at 1#for different infiltration rates

圖8 不同入滲速度時1#處溶質的物質的量濃度變化曲線Fig.8 Variations of concentration at 1#for different infiltration rates

由圖7觀察得到，土壤溫度隨時間變化大致逐漸升高，且入滲速度越大土壤溫升越快;當v=1×10－6m·s－1時，土壤溫度波動明顯.這是因為滲流初期，入滲速度越小土壤的飽和時間越長，導致上、下層土壤之間的熱導率差異縮減緩慢，從而使得土壤溫度產生波動.由圖8可觀察得到，對于同一個觀測點1#，在土壤滲流初期，滲流速度越小1#處溶質的c(A)越大;滲流后期，入滲速度越大溶質的c(A)越大.這主要是由于滲流初期土壤非飽和，溶質的對流作用大于彌散作用，當入滲速度加快時溶質對流作用越明顯，所以溶質遷移越快，溶質的c(A)增加緩慢;而滲流后期土壤達到飽和，此時彌散作用起主導作用，入滲速度越大水動力彌散系數(shù)越大，溶質的c(A)增加越快.綜合圖7，8可以得到，入滲速度的增加可以促進土壤中熱量的擴散和對流，增強溶質的對流作用和彌散作用.

3 結論

1)滲流初期，上、下層土壤含水率不同導致土壤熱導率不均勻，對土壤溫度變化具有一定影響.

2)越接近土壤地表，溶質的物質的量濃度變化受滲流影響越顯著，且隨著深度增加，影響越來越小.

3)入滲速度的增加可以促進土壤中熱量的擴散和對流，增強溶質的對流作用和彌散作用.

［1］ VIOTTIP，PAPINIM P，STRACQUALURSIN，et al.Contaminant transport in an unsaturated soil:laboratory tests and numerical simulationmodel as procedure for parameters evaluation ［J］.Ecol Model，2005，182(2):131-148.

［2］ KUMAR R P，DODAGOUDAR G R.Meshfree analysis of two-dimensional contaminant transport through unsaturated porousmedia using EFGM［J］.Int JNumer Meth Biomed Engng，2010，26(12):1797-1816.

［3］ JAVADIA A，AL-NAJJAR M M，EVANSB.Numericalmodeling of contaminant transport through soils:case study［J］.JGeotech Geoenviron Eng，2008，134(2):214-230.

［4］ GHASEMZADEH H.Heat and contaminant transport in unsaturated soil［J］.Int JCiv Eng，2008，6(2):90-107.

［5］ WINIARSKIT，LASSABATERE L，ANGULO-JARAMILLO R，et al.Characterization of the heterogeneous flow and pollutant transfer in the unsaturated zone in the fluvio-glacial deposit［J］.Procedia Environ Sci，2013，19:955-964.

［6］ 眭素剛，徐世光，劉文連，等.某渣場污染物在非飽和巖土介質中的遷移模擬研究［J］.中南大學學報(自然科學版)，2013，44(5):2173-2180.

［7］ 左自波，張璐璐，王建華.降雨條件下非飽和土中污染物遷移的數(shù)值模擬［J］.地下空間與工程學報，2011，7(S1):1347-1352.

［8］仵彥卿.多孔介質滲流與污染物遷移數(shù)學模型［M］.北京:科學出版社，2011:188-339.

［9］王華軍，路俊超，楊賓，等.溫度梯度下土壤鹽分遷移過程的數(shù)值模擬［J］.太陽能學報，2014，35(6):1086-1091.

［10］ van GENUCHTEN M T.A closed form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Sci Soc Am J，1980，44(5):892-898.

［11］ MUALEM Y.A new model for predicting the hydraulic conductivity of unsaturated porousmedia［J］.Water Resour Res，1976，12(3):513-522.

Numerical simulation of solute transport in unsaturated soils under tem perature gradient

SUN Lulu，YANGWeibo*

(Sch of Hydr，Energy ＆ Power Engin，Yangzhou Univ，Yangzhou 225127，China)

A model for solute transport in unsaturated soils is established based on the coupled governing equations of non-isothermal flow，heat transfer and solute transport.Numerical analyses are performed to investigate the effects of infiltration rates on the soil temperature variations and solute transport rules under temperature gradient.The results prove that soil temperature and solute concentration of different depths are mainly increasing over time，and both gradually reduce along the depth direction.And furthermore，the increase of the infiltration rate could promote the convection diffusion of heat in the soil and strengthen the advection and dispersion effect of solute.

unsaturated soils;temperature gradient;seepage;solute transport;numerical simulation

TK 124

A

1007-824X(2015)03-0046-05

2015-03-27.* 聯(lián)系人，E-mail:yangwb2004@163.com.

中國科學院可再生能源重點實驗室資助項目(y507k51001);江蘇省自然科學基金資助項目(BK20141278);熱流科學與工程教育部重點實驗室(西安交通大學)開放基金資助項目(KLTFSE2014KF05);廣西建筑新能源與節(jié)能重點實驗室開放基金資助項目(桂林能15-J-22-3).

孫露露，楊衛(wèi)波.溫度梯度作用下非飽和土壤中溶質遷移過程的數(shù)值模擬［J］.揚州大學學報(自然科學版)，2015，18(3):46-50.

(責任編輯 賈慧鳴，青 禾)