土壤水运动数值模拟

土壤水运动数值模拟

土壤水运动是多孔介质流运动的重要形式。土壤和地下水的各种文件过程，如大气降水、地表水入渗、地表径流、采矿蒸发和植物蒸发，与土壤水运动密切相关。对土壤水运动的数值模拟研究在大气科学、土壤学、农业工程、环境工程和地下水动力学等方面均具有重要意义。本文通过有限元法建立一维土壤水的数值模拟模型,针对土壤水运动数值模拟中对含水量的时间项比较敏感的问题,利用质量集中法进行处理,并采用特殊的迭代格式以减少模拟过程中的水量平衡误差。1土壤水价值1的模型1.1土壤含水率[]以垂直向上为正,一维土壤水的运动方程(Richard方程)可表达如下∂θ∂t=∂∂z[Κ(h)(∂h∂z+1)]-S(1)∂θ∂t=∂∂z[K(h)(∂h∂z+1)]−S(1)式中:θ为土壤体积含水率[L3L-3];h为土壤水负压[L];S为根系吸水率或其它源汇项[T-1];K(h)为饱和/非饱和水力传导度[LT-1]。本文模型采用ven-Genuchten的修正模型来描述三者之间的联系。1.2有限元方程的求解用伽辽金有限元法对式(1)进行空间离散,则对于每一个计算节点∫Ω{∂θ∂t-∂∂z[Κ∂h∂z+Κ]+S}ϕndΩ=0(2)∫Ω{∂θ∂t−∂∂z[K∂h∂z+K]+S}ϕndΩ=0(2)式中:n为计算节点编号;ϕn为线性插值基函数;Ω为一维模拟空间。将式(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分别为模拟区域上、下边界的垂向位置。将模拟区域内的土壤水负压h(z,t)用下式进行近似h′(z,t)=Ν∑n=1ϕn(z)hn(t)(4)h′(z,t)=∑n=1Nϕn(z)hn(t)(4)式中:hn(t)为待求系数,其解为各离散节点处的土壤水负压值;N为区域内的计算节点总数。在式(3)中用h′(z,t)替换h(z,t),并将积分化在各计算单元上进行可得∑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)对上式的时间项用质量集中法进行处理,即:dθndt=∑e∫Ωe∂θ∂tϕndΩ/∑e∫ΩeϕndΩ(6)dθndt=∑e∫Ωe∂θ∂tϕndΩ/∑e∫ΩeϕndΩ(6)对计算区域内的N个计算节点,应用上式(5)和式(6)进行离散处理可以得到矩阵方程组[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为节点边界流量项,Qn=Κ(∂h′∂z+1)ϕnΖΤopΖBot,对于非边界节点,由于插值基函数ϕn在上下边界处的值为0,因此Qn=0;n=m时δnm=1,n≠m时δnm=0。对于时间的离散,为保持计算的稳定性,模型中采用隐式差分的格式,即:[F]{θ}j+1-{θ}jΔtj+[A]j+1{h}j+1={Q}j+1-{B}j+1-{D}j+1(8)式中:j+1代表当前的时间层;j代表前一时间层;Δtj=tj+1-tj。模型的求解通过h进行,由于方程系数θ,A,B,D都是h的函数,因此以上方程组具有高度的非线性,在每个Δt计算时段都需要通过迭代法进行求解。迭代过程对如何处理上述矩阵方程中的含水率项十分敏感。模型中采用了一种称为“质量守恒”的方法,以减少计算过程中的水量平衡误差。这种方法在迭代过程中将式(8)中的第一项分成两部分,即:[F]{θ}j+1-{θ}jΔtj=[F]{θ}k+1j+1-{θ}kj+1Δtj+[F]{θ}kj+1-{θ}jΔtj(9)式中:k+1表示当前迭代;k表示上一次迭代。式(9)中右端第二项对于当前迭代过程来说是已知的,将式(9)中右端第一项转化为用水头表示得[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]为对角矩阵,其矩阵元素为Cnm=δnmCn,Cn为计算节点处的容水度。当迭代结果满足精度要求时,式(10)中右端第一项基本上等于0,这种特性可以有效减少求解过程中的水量平衡误差。经过上述处理之后,暂不考虑边界条件的处理,迭代求解过程中矩阵方程可以表达如下([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)当前第k+1次迭代时,式(11)右端及左端的各项系数矩阵均已知,经整理可得编程求解形式[G]{h}={g}(12)式中:G为系数矩阵;{g}为右端已知项。如果节点编号由上自下顺序进行,则G为主对角线占优的三对角矩阵,采用追赶法求解。1.3地表蒸发状态土壤水方程的定解条件有两个:一是初始条件h(z,0)=h0(z),二是边界条件。即使是在一维情况下,土壤水的上边界条件也具有相当的复杂性。土壤水的上边界主要接受降雨/灌溉,以及土表蒸发,受限于土表的入渗能力,在降雨/灌溉过程中由于强度的大小不同,土表可能积水或不积水,也可能超出积水滞蓄深度产流。在地表处于蒸发状态时,根据土壤输水能力,也可分为达到或未达到蒸发极限的情况。目前大多数土壤水模型对于土表积水(SurfacePonding)过程以及产流过程尚无严格的数值方法,通常积水和非积水过程人为分离模拟,本文试图从边界处理的角度解决这一问题。1.3.1土壤水负压h当灌溉/降雨强度,或蒸发强度并未超出土壤入渗/蒸发能力时,上边界条件为流量边界条件(二类边界条件)-Κ(∂h/∂z+1)=Ι(t)‚hA<h<0(13)式中:hA为对应土壤风干含水率时的土壤水负压,由土质而定;I(t)为作用在土表之上的潜在净流入/流出强度。由于模型取z轴向上为正,I(t)在数值上等于t时刻的外界蒸发强度减去降雨/灌溉强度。二类边界条件可以直接应用,在上边界节点处,有Q=Κ(∂h′∂z+1)ϕΖΤopΖΤopΖBot=-ˉΙΔtj,其中ϕZTop为上边界节点处的线性插值基函数,ˉΙΔtj为计算时段Δtj内的平均潜在净流入/流出强度(下同)。作为己知项,-ˉΙΔtj可以并入式(12)的右端项[g]中。1.3.2tj边界节点的流量强度计算当降雨/灌溉强度超出土壤入渗能力时,地表将会出现积水现象,积水深度在没有超出最大深度(对于田间灌溉可以理解为田埂高度,对于降雨产流则可以理解为地表洼地滞蓄深度)之前,超渗的水量将会在地表随时间积累。对于一维问题,该边界条件可用下式描述-∂h/∂t=Ι(t)+Κ(∂h/∂z+1)‚0≤h≤hS(14)式中:hS为地表最大积水滞蓄深度。这种边界情况下h既是地表的负压水头(此时为静水压力),又代表了地表的积水深度。由于此时边界表达式中含有时间项,该边界形式不同于一般的流量或水头边界,本文称之为积水边界条件。应用该边界条件时需要对时间进行离散。为保持与土壤水控制方程一致,离散时同样采用隐式格式。对于计算时段Δtj,有Q=Κ(∂h∂z+1)ϕΖΤop|ΖΤopΖBot=-[Ι(t)+∂h∂t)ϕΖΤop]|ΖΤopΖBot=-ˉΙΔtj-hj+1ΖΤop-hjΖΤopΔtj(15)上式中含有未知项hj+1ΖΤop,即上边界节点处当前待求的水头值,其系数1/Δtj可以并入式(12)的系数矩阵[G]的主对角线中,而已知项-(ˉΙΔtj+hjΖΤop/Δtj)则可以并入式(12)的右端项{g}中。式(15)具有较明显的物理意义,右端第一项为Δtj时间内上边界作用的平均潜在净流入/流出强度,第二项为Δtj时间内上边界积水深度从hjΖΤop变化到hj+1ΖΤop引起的积水量变化强度,两者的综合叠加结果即为上边界通过的实际流量强度。当Δtj时段计算收敛后,即已经求出tj+1时刻各计算节点处的土壤水负压{h}j+1,则可以通过上边界处的节点方程确定上边界的实际流量强度Q。如果在模型中n代表上边界节点的编号,由式(8)可知Δtj时段内上边界节点处的实际流量强度为Qn=Ν∑m=1Fnmθj+1m-θjmΔtj+Ν∑m=1Aj+1nmhj+1m+Bj+1n+Dj+1n(16)其中:N为模拟区域内的计算节点个数。式(15)适用于地表已经积水时的情况,即在tj和tj+1时刻,上边界节点处的土壤水负压hZTop均大于0(静水压力)。但在边界积水过程中有两个特殊的时段需要单独处理,一是上边界从非积水状态过渡到积水状态时(hjΖΤop<0,hj+1ΖΤop≥0),二是从积水状态过渡到非积水状态时(hjΖΤop≥0,hj+1ΖΤop<0)。从边界上的水量平衡出发,结合式(15)的物理意义。对于第一种情况,在该Δtj时段内上边界实际的流量强度为Q=-ˉΙΔtj-(hj+1ΖΤop-0)/Δtj(17)上式右端第二项为Δtj时间内边界上的积水厚度从0变化到hj+1ΖΤop引起的积水量变化强度。与积水状态时的处理相同,系数1Δtj并入式(12)的系数矩阵[G],而已知项-ˉΙΔtj则可以并入式(12)的右端项{g}。在该计算时段求解完毕后,边界上实际通过的流量大小仍可由式(16)计算。对于第二种情况,从上边界的水量平衡观点出发,此时边界上的实际流量强度可以表达为Q=-ˉΙΔtj-(0-hjΖΤop)/Δtj(18)上式右端第二项为Δtj时间内边界上的积水厚度从hjΖΤop变化到0引起的积水量变化强度。注意到对于当前时段Δtj,边界节点在tj时刻的积水深度hjΖΤop已知,因此在该时段内Q为已知的,将-ˉΙΔtj+hjΖΤop/Δtj并入式(12)中的右端项{g}即可。1.3.3建立边界模型当积水深度超出地表最大滞蓄深度hS时,积水深度不再增加,超渗的水量即形成地表径流。此时的上边界条件可以用水头边界条件(一类边界)描述。h=hS,h>hS(19)为了模拟上边界积水深度达到极限深度后不再增加的情况,在数值模拟过程中可如下处理:计算时段Δtj内每次应用上述的积水边界条件进行计算时,迭代计算收敛之后都判断一次积水深度hj+1ΖΤop的大小,如果超过极限深度hS,则将该时段的计算结果作废,上边界处的节点方程用h=hS代替,并对该时段进行重新计算。该时段计算完毕后上边界的实际入渗强度Qn可由式(16)确定。根据上边界处的水量平衡关系,可得Δt时段内土表的产流量为RΔtj=-[-ˉΙΔtj-Qn]Δt。1.3.4土壤水负压h在外界的蒸发力超过土壤的输水能力时,土表处于风干状态,土表的实际蒸发强度与外界蒸发力的大小将不再存在直接关系,此时的边界条件为水头边界,可描述为h=hA‚h≥hA(20)式中:hA为对应土壤风干含水率时的土壤水负压。在满足式(20)的判断条件时,迭代过程中上边界处的节点方程将由h=hA代替,上边界的实际蒸发强度可由式(16)确定。2模型验证2.1土壤水运动特性分析模型的正确性首先用试验室一维入渗试验验证,二维有限元土壤水模型SWMS—2D曾利用该入渗试验作过模型模拟测试。本文模型对该试验进行同样的数值模拟,并将模拟结果与SWMS—2D的模拟结果进行对比。该入渗试验过程如下(图1):一维砂柱的高度为61cm,初始土壤水负压为-150cm;砂柱中的土质均一,各向同性,饱和渗透系数为0.043cm/min;试验过程中土柱顶端维持2cm厚的水层,底端封闭不排水;试验持续的时间为90min。本次数值模拟时的计算节点总数为56个,计算单元55个,节点编号和单元编号的分布情况见图1,所取土壤水运动参数见表1(与SWMS—2D所取参数一致)。模拟过程中上边界应用定水头边界(h=2cm),下边界应用隔水边界(Q=0)。数值模拟过程中土柱剖面土壤含水率变化情况如图2所示。从本模型90min时模拟的土壤含水率分布与SWMS—2D模拟的结果对比来看,两模型之间模拟的结果十分接近,说明本文模型在基本计算方面是正确的。2.2砂柱剖面含水率结果分析以上用简单的一维入渗模拟验证了模型的基本正确性,为了验证在降雨/灌溉、蒸发强度给定的情况下,上边界对积水过程的模拟能力,维持以上模拟环境中的土壤水运动参数、初始条件和单元剖分情况不变,仅对上、下边界的情况进行修改。模拟的总时长为480min,对于上边界,假设在0～60min时间内作用0.5cm/min的降雨强度,60～240min时间内无降雨、蒸发,240～480min时间内作用0.0625cm/min的蒸发强度;对于下边界,假设0～240min时间内封闭不排水(Q=0),240～480min时间内自由排水(排水期间下边界节点处h=0)。实际情况下不可能出现如此大强度的降雨和蒸发,这里只是作为检验模型所用。模拟过程中暂不限制上边界积水深度(hS取很大的值)。首先观察上、下边界节点在模拟过程中的变化情况。由图3可知,在0～60min的降雨过程中,由于降雨强度很大,上边界节点在计算过程中迅速饱和并开始积水,在降雨结束时刻(60min)边界积水深度达到最大值21.27cm。随后60～240min之内上边界无降雨、蒸发,且下边界无排水,然而在160min之前积水深度减少,这是因为砂柱内部尚未完全饱和,积水可继续入渗。160min之后砂柱完全饱和,积水深度维持13.33cm不变直到第240min。240～480min之内上边界开始有蒸发作用,同时下边界开始排水,因此上边界积水深度开始下降。在360min以后上边界已无积水,节点处的土壤水负压在高强度的蒸发条件下迅速增加。由图4可知,由于在240min之前底边界隔水,在上边界的入渗影响尚未到达砂柱底端之前,下边界节点处的土壤水负压变化很小,115min之后入渗水量到达下边界,该边界节点的负压迅速减小并向正值过渡,160min时砂柱完全饱和,下边界节点处的土壤水负压达到正的最大值74.33cm,并维持该值一直到240min之后下边界开始自由排水,此时模型将下边界切换为定水头边界(h=0)并维持到模拟结束。模拟过程中不同时刻的砂柱剖面含水率分布见图5。从图中分析可知边界水头的变化与剖面含水率之间存在良好对应关系。模拟过程中上、下边界处的累计流入/流出水量见图6,以流入为正,流出为负。由图5可知,对于上边界,在160min之前上边界持续有水量入渗,水量持续正累积。160min之后砂柱完全饱和,上边界虽然有积水(图3)但水量不能入渗,因此在160～240min之间累积水量维持16.7cm3无变化。240min之后下边界开始排水,上边界的积水又开始入渗。在第360min时积水在入渗和上边界蒸发的双重作用下消耗完毕,此时上边界已经正累积入渗22.4cm3的水量,此后上边界在蒸发作用之下开始流量负累积直到模拟结束。对于下边界,由于0～240min之内为隔水边界,因此累积流入/流出水量维持为0;240min之后边界开始自由排水,边界上持续有水量负累积。从以上模拟结果分析中可知,模拟过程中边界处的土壤水负压变化曲线、边界的流入/流出水量累积曲线及各时刻砂柱剖面含水率分布之间均存在良好对应关系,从一方面证明了模型的正确性,然而模拟过程中各个时刻的水量是否平衡、有无明显计算误差尚需分析确定。根据水量平衡关系,对于本次模拟应该有:(1)在上边界积水过程中任意时刻(0～360min):累计降雨量=累计蒸发量+上边界累计入渗水量+上边界积水量;(2)在数值模拟过程中任意时刻(0～480min):砂柱含水量=砂柱初始含水量+上边界累积流入/流出水量+下边界累积流入/流出水量。在模拟过程中选取不同的时刻进行分析,其结果见表2和表3。分析结果表明模拟计算过程中以上两种水量平衡关系确实成立,而且计算误差很小。以上积水过程验证并未考虑积水滞蓄深度限