土壤水運(yùn)動(dòng)數(shù)值模擬

土壤水運(yùn)動(dòng)數(shù)值模擬

土壤水運(yùn)動(dòng)是多孔介質(zhì)流運(yùn)動(dòng)的重要形式。土壤和地下水的各種文件過(guò)程，如大氣降水、地表水入滲、地表徑流、采礦蒸發(fā)和植物蒸發(fā)，與土壤水運(yùn)動(dòng)密切相關(guān)。對(duì)土壤水運(yùn)動(dòng)的數(shù)值模擬研究在大氣科學(xué)、土壤學(xué)、農(nóng)業(yè)工程、環(huán)境工程和地下水動(dòng)力學(xué)等方面均具有重要意義。本文通過(guò)有限元法建立一維土壤水的數(shù)值模擬模型,針對(duì)土壤水運(yùn)動(dòng)數(shù)值模擬中對(duì)含水量的時(shí)間項(xiàng)比較敏感的問(wèn)題,利用質(zhì)量集中法進(jìn)行處理,并采用特殊的迭代格式以減少模擬過(guò)程中的水量平衡誤差。1土壤水價(jià)值1的模型1.1土壤含水率[]以垂直向上為正,一維土壤水的運(yùn)動(dòng)方程(Richard方程)可表達(dá)如下?θ?t=??z[Κ(h)(?h?z+1)]-S(1)?θ?t=??z[K(h)(?h?z+1)]?S(1)式中:θ為土壤體積含水率[L3L-3];h為土壤水負(fù)壓[L];S為根系吸水率或其它源匯項(xiàng)[T-1];K(h)為飽和/非飽和水力傳導(dǎo)度[LT-1]。本文模型采用ven-Genuchten的修正模型來(lái)描述三者之間的聯(lián)系。1.2有限元方程的求解用伽遼金有限元法對(duì)式(1)進(jìn)行空間離散,則對(duì)于每一個(gè)計(jì)算節(jié)點(diǎn)∫Ω{?θ?t-??z[Κ?h?z+Κ]+S}?ndΩ=0(2)∫Ω{?θ?t???z[K?h?z+K]+S}?ndΩ=0(2)式中:n為計(jì)算節(jié)點(diǎn)編號(hào);?n為線性插值基函數(shù);Ω為一維模擬空間。將式(2)展開并利用分部積分法,有∫Ω?θ?t?ndΩ+∫Ω(Κ?h?z+Κ)??n?zdΩ-(Κ?h?z+Κ)?n|ΖΤopΖBot+∫ΩS?ndΩ=0(3)∫Ω?θ?t?ndΩ+∫Ω(K?h?z+K)??n?zdΩ?(K?h?z+K)?n|ZTopZBot+∫ΩS?ndΩ=0(3)式中:ZTop和ZBot分別為模擬區(qū)域上、下邊界的垂向位置。將模擬區(qū)域內(nèi)的土壤水負(fù)壓h(z,t)用下式進(jìn)行近似h′(z,t)=Ν∑n=1?n(z)hn(t)(4)h′(z,t)=∑n=1N?n(z)hn(t)(4)式中:hn(t)為待求系數(shù),其解為各離散節(jié)點(diǎn)處的土壤水負(fù)壓值;N為區(qū)域內(nèi)的計(jì)算節(jié)點(diǎn)總數(shù)。在式(3)中用h′(z,t)替換h(z,t),并將積分化在各計(jì)算單元上進(jìn)行可得∑e∫Ωe?θ?t?ndΩ+∑e∫ΩeΚ?h′?z??n?zdΩ=Κ(?h′?z+1)?n|ΖΤopΖBot-∑e∫ΩeΚ??n?zdΩ-∑e∫ΩeS?ndΩ(5)∑e∫Ωe?θ?t?ndΩ+∑e∫ΩeK?h′?z??n?zdΩ=K(?h′?z+1)?n|ZTopZBot?∑e∫ΩeK??n?zdΩ?∑e∫ΩeS?ndΩ(5)對(duì)上式的時(shí)間項(xiàng)用質(zhì)量集中法進(jìn)行處理,即:dθndt=∑e∫Ωe?θ?t?ndΩ/∑e∫Ωe?ndΩ(6)dθndt=∑e∫Ωe?θ?t?ndΩ/∑e∫Ωe?ndΩ(6)對(duì)計(jì)算區(qū)域內(nèi)的N個(gè)計(jì)算節(jié)點(diǎn),應(yīng)用上式(5)和式(6)進(jìn)行離散處理可以得到矩陣方程組[F]d{θ}dt+[A]{h}={Q}-{B}-{D}(7)式中:Fnm=δnm∑e∫Ωe?ndΩ;Anm=∑e∫ΩeΚ??m?z??n?zdΩ;Bn=∑e∫ΩeΚ??n?zdΩ;Dn=∑e∫ΩeS?ndΩ;Ωn為節(jié)點(diǎn)邊界流量項(xiàng),Qn=Κ(?h′?z+1)?nΖΤopΖBot,對(duì)于非邊界節(jié)點(diǎn),由于插值基函數(shù)?n在上下邊界處的值為0,因此Qn=0;n=m時(shí)δnm=1,n≠m時(shí)δnm=0。對(duì)于時(shí)間的離散,為保持計(jì)算的穩(wěn)定性,模型中采用隱式差分的格式,即:[F]{θ}j+1-{θ}jΔtj+[A]j+1{h}j+1={Q}j+1-{B}j+1-{D}j+1(8)式中:j+1代表當(dāng)前的時(shí)間層;j代表前一時(shí)間層;Δtj=tj+1-tj。模型的求解通過(guò)h進(jìn)行,由于方程系數(shù)θ,A,B,D都是h的函數(shù),因此以上方程組具有高度的非線性,在每個(gè)Δt計(jì)算時(shí)段都需要通過(guò)迭代法進(jìn)行求解。迭代過(guò)程對(duì)如何處理上述矩陣方程中的含水率項(xiàng)十分敏感。模型中采用了一種稱為“質(zhì)量守恒”的方法,以減少計(jì)算過(guò)程中的水量平衡誤差。這種方法在迭代過(guò)程中將式(8)中的第一項(xiàng)分成兩部分,即:[F]{θ}j+1-{θ}jΔtj=[F]{θ}k+1j+1-{θ}kj+1Δtj+[F]{θ}kj+1-{θ}jΔtj(9)式中:k+1表示當(dāng)前迭代;k表示上一次迭代。式(9)中右端第二項(xiàng)對(duì)于當(dāng)前迭代過(guò)程來(lái)說(shuō)是已知的,將式(9)中右端第一項(xiàng)轉(zhuǎn)化為用水頭表示得[F]{θ}j+1-{θ}jΔtj=[F][C]kj+1{h}k+1j+1-{h}kj+1Δtj+[F]{θ}kj+1-{θ}jΔtj(10)式中:[C]、[F]為對(duì)角矩陣,其矩陣元素為Cnm=δnmCn,Cn為計(jì)算節(jié)點(diǎn)處的容水度。當(dāng)?shù)Y(jié)果滿足精度要求時(shí),式(10)中右端第一項(xiàng)基本上等于0,這種特性可以有效減少求解過(guò)程中的水量平衡誤差。經(jīng)過(guò)上述處理之后,暫不考慮邊界條件的處理,迭代求解過(guò)程中矩陣方程可以表達(dá)如下([F][C]kj+1Δtj+[A]kj+1){h}k+1j+1=[F][C]kj+1Δtj{h}kj+1-[F]{θ}kj+1-{θ}jΔtj-{B}kj+1-{D}kj+1(11)當(dāng)前第k+1次迭代時(shí),式(11)右端及左端的各項(xiàng)系數(shù)矩陣均已知,經(jīng)整理可得編程求解形式[G]{h}={g}(12)式中:G為系數(shù)矩陣;{g}為右端已知項(xiàng)。如果節(jié)點(diǎn)編號(hào)由上自下順序進(jìn)行,則G為主對(duì)角線占優(yōu)的三對(duì)角矩陣,采用追趕法求解。1.3地表蒸發(fā)狀態(tài)土壤水方程的定解條件有兩個(gè):一是初始條件h(z,0)=h0(z),二是邊界條件。即使是在一維情況下,土壤水的上邊界條件也具有相當(dāng)?shù)膹?fù)雜性。土壤水的上邊界主要接受降雨/灌溉,以及土表蒸發(fā),受限于土表的入滲能力,在降雨/灌溉過(guò)程中由于強(qiáng)度的大小不同,土表可能積水或不積水,也可能超出積水滯蓄深度產(chǎn)流。在地表處于蒸發(fā)狀態(tài)時(shí),根據(jù)土壤輸水能力,也可分為達(dá)到或未達(dá)到蒸發(fā)極限的情況。目前大多數(shù)土壤水模型對(duì)于土表積水(SurfacePonding)過(guò)程以及產(chǎn)流過(guò)程尚無(wú)嚴(yán)格的數(shù)值方法,通常積水和非積水過(guò)程人為分離模擬,本文試圖從邊界處理的角度解決這一問(wèn)題。1.3.1土壤水負(fù)壓h當(dāng)灌溉/降雨強(qiáng)度,或蒸發(fā)強(qiáng)度并未超出土壤入滲/蒸發(fā)能力時(shí),上邊界條件為流量邊界條件(二類邊界條件)-Κ(?h/?z+1)=Ι(t)?hA<h<0(13)式中:hA為對(duì)應(yīng)土壤風(fēng)干含水率時(shí)的土壤水負(fù)壓,由土質(zhì)而定;I(t)為作用在土表之上的潛在凈流入/流出強(qiáng)度。由于模型取z軸向上為正,I(t)在數(shù)值上等于t時(shí)刻的外界蒸發(fā)強(qiáng)度減去降雨/灌溉強(qiáng)度。二類邊界條件可以直接應(yīng)用,在上邊界節(jié)點(diǎn)處,有Q=Κ(?h′?z+1)?ΖΤopΖΤopΖBot=-ˉΙΔtj,其中?ZTop為上邊界節(jié)點(diǎn)處的線性插值基函數(shù),ˉΙΔtj為計(jì)算時(shí)段Δtj內(nèi)的平均潛在凈流入/流出強(qiáng)度(下同)。作為己知項(xiàng),-ˉΙΔtj可以并入式(12)的右端項(xiàng)[g]中。1.3.2tj邊界節(jié)點(diǎn)的流量強(qiáng)度計(jì)算當(dāng)降雨/灌溉強(qiáng)度超出土壤入滲能力時(shí),地表將會(huì)出現(xiàn)積水現(xiàn)象,積水深度在沒(méi)有超出最大深度(對(duì)于田間灌溉可以理解為田埂高度,對(duì)于降雨產(chǎn)流則可以理解為地表洼地滯蓄深度)之前,超滲的水量將會(huì)在地表隨時(shí)間積累。對(duì)于一維問(wèn)題,該邊界條件可用下式描述-?h/?t=Ι(t)+Κ(?h/?z+1)?0≤h≤hS(14)式中:hS為地表最大積水滯蓄深度。這種邊界情況下h既是地表的負(fù)壓水頭(此時(shí)為靜水壓力),又代表了地表的積水深度。由于此時(shí)邊界表達(dá)式中含有時(shí)間項(xiàng),該邊界形式不同于一般的流量或水頭邊界,本文稱之為積水邊界條件。應(yīng)用該邊界條件時(shí)需要對(duì)時(shí)間進(jìn)行離散。為保持與土壤水控制方程一致,離散時(shí)同樣采用隱式格式。對(duì)于計(jì)算時(shí)段Δtj,有Q=Κ(?h?z+1)?ΖΤop|ΖΤopΖBot=-[Ι(t)+?h?t)?ΖΤop]|ΖΤopΖBot=-ˉΙΔtj-hj+1ΖΤop-hjΖΤopΔtj(15)上式中含有未知項(xiàng)hj+1ΖΤop,即上邊界節(jié)點(diǎn)處當(dāng)前待求的水頭值,其系數(shù)1/Δtj可以并入式(12)的系數(shù)矩陣[G]的主對(duì)角線中,而已知項(xiàng)-(ˉΙΔtj+hjΖΤop/Δtj)則可以并入式(12)的右端項(xiàng){g}中。式(15)具有較明顯的物理意義,右端第一項(xiàng)為Δtj時(shí)間內(nèi)上邊界作用的平均潛在凈流入/流出強(qiáng)度,第二項(xiàng)為Δtj時(shí)間內(nèi)上邊界積水深度從hjΖΤop變化到hj+1ΖΤop引起的積水量變化強(qiáng)度,兩者的綜合疊加結(jié)果即為上邊界通過(guò)的實(shí)際流量強(qiáng)度。當(dāng)Δtj時(shí)段計(jì)算收斂后,即已經(jīng)求出tj+1時(shí)刻各計(jì)算節(jié)點(diǎn)處的土壤水負(fù)壓{h}j+1,則可以通過(guò)上邊界處的節(jié)點(diǎn)方程確定上邊界的實(shí)際流量強(qiáng)度Q。如果在模型中n代表上邊界節(jié)點(diǎn)的編號(hào),由式(8)可知Δtj時(shí)段內(nèi)上邊界節(jié)點(diǎn)處的實(shí)際流量強(qiáng)度為Qn=Ν∑m=1Fnmθj+1m-θjmΔtj+Ν∑m=1Aj+1nmhj+1m+Bj+1n+Dj+1n(16)其中:N為模擬區(qū)域內(nèi)的計(jì)算節(jié)點(diǎn)個(gè)數(shù)。式(15)適用于地表已經(jīng)積水時(shí)的情況,即在tj和tj+1時(shí)刻,上邊界節(jié)點(diǎn)處的土壤水負(fù)壓hZTop均大于0(靜水壓力)。但在邊界積水過(guò)程中有兩個(gè)特殊的時(shí)段需要單獨(dú)處理,一是上邊界從非積水狀態(tài)過(guò)渡到積水狀態(tài)時(shí)(hjΖΤop<0,hj+1ΖΤop≥0),二是從積水狀態(tài)過(guò)渡到非積水狀態(tài)時(shí)(hjΖΤop≥0,hj+1ΖΤop<0)。從邊界上的水量平衡出發(fā),結(jié)合式(15)的物理意義。對(duì)于第一種情況,在該Δtj時(shí)段內(nèi)上邊界實(shí)際的流量強(qiáng)度為Q=-ˉΙΔtj-(hj+1ΖΤop-0)/Δtj(17)上式右端第二項(xiàng)為Δtj時(shí)間內(nèi)邊界上的積水厚度從0變化到hj+1ΖΤop引起的積水量變化強(qiáng)度。與積水狀態(tài)時(shí)的處理相同,系數(shù)1Δtj并入式(12)的系數(shù)矩陣[G],而已知項(xiàng)-ˉΙΔtj則可以并入式(12)的右端項(xiàng){g}。在該計(jì)算時(shí)段求解完畢后,邊界上實(shí)際通過(guò)的流量大小仍可由式(16)計(jì)算。對(duì)于第二種情況,從上邊界的水量平衡觀點(diǎn)出發(fā),此時(shí)邊界上的實(shí)際流量強(qiáng)度可以表達(dá)為Q=-ˉΙΔtj-(0-hjΖΤop)/Δtj(18)上式右端第二項(xiàng)為Δtj時(shí)間內(nèi)邊界上的積水厚度從hjΖΤop變化到0引起的積水量變化強(qiáng)度。注意到對(duì)于當(dāng)前時(shí)段Δtj,邊界節(jié)點(diǎn)在tj時(shí)刻的積水深度hjΖΤop已知,因此在該時(shí)段內(nèi)Q為已知的,將-ˉΙΔtj+hjΖΤop/Δtj并入式(12)中的右端項(xiàng){g}即可。1.3.3建立邊界模型當(dāng)積水深度超出地表最大滯蓄深度hS時(shí),積水深度不再增加,超滲的水量即形成地表徑流。此時(shí)的上邊界條件可以用水頭邊界條件(一類邊界)描述。h=hS,h>hS(19)為了模擬上邊界積水深度達(dá)到極限深度后不再增加的情況,在數(shù)值模擬過(guò)程中可如下處理:計(jì)算時(shí)段Δtj內(nèi)每次應(yīng)用上述的積水邊界條件進(jìn)行計(jì)算時(shí),迭代計(jì)算收斂之后都判斷一次積水深度hj+1ΖΤop的大小,如果超過(guò)極限深度hS,則將該時(shí)段的計(jì)算結(jié)果作廢,上邊界處的節(jié)點(diǎn)方程用h=hS代替,并對(duì)該時(shí)段進(jìn)行重新計(jì)算。該時(shí)段計(jì)算完畢后上邊界的實(shí)際入滲強(qiáng)度Qn可由式(16)確定。根據(jù)上邊界處的水量平衡關(guān)系,可得Δt時(shí)段內(nèi)土表的產(chǎn)流量為RΔtj=-[-ˉΙΔtj-Qn]Δt。1.3.4土壤水負(fù)壓h在外界的蒸發(fā)力超過(guò)土壤的輸水能力時(shí),土表處于風(fēng)干狀態(tài),土表的實(shí)際蒸發(fā)強(qiáng)度與外界蒸發(fā)力的大小將不再存在直接關(guān)系,此時(shí)的邊界條件為水頭邊界,可描述為h=hA?h≥hA(20)式中:hA為對(duì)應(yīng)土壤風(fēng)干含水率時(shí)的土壤水負(fù)壓。在滿足式(20)的判斷條件時(shí),迭代過(guò)程中上邊界處的節(jié)點(diǎn)方程將由h=hA代替,上邊界的實(shí)際蒸發(fā)強(qiáng)度可由式(16)確定。2模型驗(yàn)證2.1土壤水運(yùn)動(dòng)特性分析模型的正確性首先用試驗(yàn)室一維入滲試驗(yàn)驗(yàn)證,二維有限元土壤水模型SWMS—2D曾利用該入滲試驗(yàn)作過(guò)模型模擬測(cè)試。本文模型對(duì)該試驗(yàn)進(jìn)行同樣的數(shù)值模擬,并將模擬結(jié)果與SWMS—2D的模擬結(jié)果進(jìn)行對(duì)比。該入滲試驗(yàn)過(guò)程如下(圖1):一維砂柱的高度為61cm,初始土壤水負(fù)壓為-150cm;砂柱中的土質(zhì)均一,各向同性,飽和滲透系數(shù)為0.043cm/min;試驗(yàn)過(guò)程中土柱頂端維持2cm厚的水層,底端封閉不排水;試驗(yàn)持續(xù)的時(shí)間為90min。本次數(shù)值模擬時(shí)的計(jì)算節(jié)點(diǎn)總數(shù)為56個(gè),計(jì)算單元55個(gè),節(jié)點(diǎn)編號(hào)和單元編號(hào)的分布情況見圖1,所取土壤水運(yùn)動(dòng)參數(shù)見表1(與SWMS—2D所取參數(shù)一致)。模擬過(guò)程中上邊界應(yīng)用定水頭邊界(h=2cm),下邊界應(yīng)用隔水邊界(Q=0)。數(shù)值模擬過(guò)程中土柱剖面土壤含水率變化情況如圖2所示。從本模型90min時(shí)模擬的土壤含水率分布與SWMS—2D模擬的結(jié)果對(duì)比來(lái)看,兩模型之間模擬的結(jié)果十分接近,說(shuō)明本文模型在基本計(jì)算方面是正確的。2.2砂柱剖面含水率結(jié)果分析以上用簡(jiǎn)單的一維入滲模擬驗(yàn)證了模型的基本正確性,為了驗(yàn)證在降雨/灌溉、蒸發(fā)強(qiáng)度給定的情況下,上邊界對(duì)積水過(guò)程的模擬能力,維持以上模擬環(huán)境中的土壤水運(yùn)動(dòng)參數(shù)、初始條件和單元剖分情況不變,僅對(duì)上、下邊界的情況進(jìn)行修改。模擬的總時(shí)長(zhǎng)為480min,對(duì)于上邊界,假設(shè)在0～60min時(shí)間內(nèi)作用0.5cm/min的降雨強(qiáng)度,60～240min時(shí)間內(nèi)無(wú)降雨、蒸發(fā),240～480min時(shí)間內(nèi)作用0.0625cm/min的蒸發(fā)強(qiáng)度;對(duì)于下邊界,假設(shè)0～240min時(shí)間內(nèi)封閉不排水(Q=0),240～480min時(shí)間內(nèi)自由排水(排水期間下邊界節(jié)點(diǎn)處h=0)。實(shí)際情況下不可能出現(xiàn)如此大強(qiáng)度的降雨和蒸發(fā),這里只是作為檢驗(yàn)?zāi)Ｐ退?。模擬過(guò)程中暫不限制上邊界積水深度(hS取很大的值)。首先觀察上、下邊界節(jié)點(diǎn)在模擬過(guò)程中的變化情況。由圖3可知,在0～60min的降雨過(guò)程中,由于降雨強(qiáng)度很大,上邊界節(jié)點(diǎn)在計(jì)算過(guò)程中迅速飽和并開始積水,在降雨結(jié)束時(shí)刻(60min)邊界積水深度達(dá)到最大值21.27cm。隨后60～240min之內(nèi)上邊界無(wú)降雨、蒸發(fā),且下邊界無(wú)排水,然而在160min之前積水深度減少,這是因?yàn)樯爸鶅?nèi)部尚未完全飽和,積水可繼續(xù)入滲。160min之后砂柱完全飽和,積水深度維持13.33cm不變直到第240min。240～480min之內(nèi)上邊界開始有蒸發(fā)作用,同時(shí)下邊界開始排水,因此上邊界積水深度開始下降。在360min以后上邊界已無(wú)積水,節(jié)點(diǎn)處的土壤水負(fù)壓在高強(qiáng)度的蒸發(fā)條件下迅速增加。由圖4可知,由于在240min之前底邊界隔水,在上邊界的入滲影響尚未到達(dá)砂柱底端之前,下邊界節(jié)點(diǎn)處的土壤水負(fù)壓變化很小,115min之后入滲水量到達(dá)下邊界,該邊界節(jié)點(diǎn)的負(fù)壓迅速減小并向正值過(guò)渡,160min時(shí)砂柱完全飽和,下邊界節(jié)點(diǎn)處的土壤水負(fù)壓達(dá)到正的最大值74.33cm,并維持該值一直到240min之后下邊界開始自由排水,此時(shí)模型將下邊界切換為定水頭邊界(h=0)并維持到模擬結(jié)束。模擬過(guò)程中不同時(shí)刻的砂柱剖面含水率分布見圖5。從圖中分析可知邊界水頭的變化與剖面含水率之間存在良好對(duì)應(yīng)關(guān)系。模擬過(guò)程中上、下邊界處的累計(jì)流入/流出水量見圖6,以流入為正,流出為負(fù)。由圖5可知,對(duì)于上邊界,在160min之前上邊界持續(xù)有水量入滲,水量持續(xù)正累積。160min之后砂柱完全飽和,上邊界雖然有積水(圖3)但水量不能入滲,因此在160～240min之間累積水量維持16.7cm3無(wú)變化。240min之后下邊界開始排水,上邊界的積水又開始入滲。在第360min時(shí)積水在入滲和上邊界蒸發(fā)的雙重作用下消耗完畢,此時(shí)上邊界已經(jīng)正累積入滲22.4cm3的水量,此后上邊界在蒸發(fā)作用之下開始流量負(fù)累積直到模擬結(jié)束。對(duì)于下邊界,由于0～240min之內(nèi)為隔水邊界,因此累積流入/流出水量維持為0;240min之后邊界開始自由排水,邊界上持續(xù)有水量負(fù)累積。從以上模擬結(jié)果分析中可知,模擬過(guò)程中邊界處的土壤水負(fù)壓變化曲線、邊界的流入/流出水量累積曲線及各時(shí)刻砂柱剖面含水率分布之間均存在良好對(duì)應(yīng)關(guān)系,從一方面證明了模型的正確性,然而模擬過(guò)程中各個(gè)時(shí)刻的水量是否平衡、有無(wú)明顯計(jì)算誤差尚需分析確定。根據(jù)水量平衡關(guān)系,對(duì)于本次模擬應(yīng)該有:(1)在上邊界積水過(guò)程中任意時(shí)刻(0～360min):累計(jì)降雨量=累計(jì)蒸發(fā)量+上邊界累計(jì)入滲水量+上邊界積水量;(2)在數(shù)值模擬過(guò)程中任意時(shí)刻(0～480min):砂柱含水量=砂柱初始含水量+上邊界累積流入/流出水量+下邊界累積流入/流出水量。在模擬過(guò)程中選取不同的時(shí)刻進(jìn)行分析,其結(jié)果見表2和表3。分析結(jié)果表明模擬計(jì)算過(guò)程中以上兩種水量平衡關(guān)系確實(shí)成立,而且計(jì)算誤差很小。以上積水過(guò)程驗(yàn)證并未考慮積水滯蓄深度限