1997年5月水 利 学 报SHUILI XUEBAO第5期地下水与土壤盐渍化关系的动态模拟陈亚新 史海滨 田存旺(内蒙古农牧学院水利工程系)提要在大气-土壤-地下水连续系统(ASGC)中,本文将土壤水-地下水部分用动力学理论和数值模拟技术结合的方法,研究了地下水与土壤盐渍化关系,并进行了动态模拟.在完成土壤水盐数据自动采集系统研制、水盐运动参数测定、数学模型建模识别、数值计算方法选择评价4项基础研究后,针对干旱区的一个大型灌区土壤-水环境中溶质迁移动态进行了系统模拟,获得复杂环境条件下多因子、多处理组合模式化的图谱,有助地下水与土壤盐渍化定量关系深入揭示、水土资源管理与预测.关键词 地下水,土壤盐渍化,动态模拟,溶质迁移,相关性.一、导 言大气-土壤-地下水连续体(Atmosphere2Soil2Groundwater Continuum简称ASGC)中,土壤盐渍化与地下水文关系密切,气候条件和矿质地下水动态是决定盐渍化敏感的环境因子、是包气带内溶质迁移能量转换和表土积盐的重要机制.地下水与土壤盐渍化关系研究过程的发展及其代表性的研究方法有以下4种: (1)五十年代,主要是应用调查统计法.通过大面积土壤和水文地质调查,统计地下水埋深(矿化度)与表土盐渍化现象的关系,获得地下水临界深度HK.但不能精确区分地下水水质,水位动态的具体影响,不涉及溶质迁移的动力机制,所称“表土盐渍化”难以量化,只能作自然现象的 宏观描述.(2)六十年代我国用土壤强烈毛管水上升高度(毛管破裂点)加主要根系活动层深度,计算HK(或安全深度) .从构成HK的微观揭示和量化方面较国外传统理论有较大进步1.8. ,仍不能对饱和2非饱和流连续体的积盐和脱盐能量转化过程进行深入揭示.(3)七十至八十年代,地中渗透仪(Lysimeters)方法的广泛应用,可进行补排关系潜水蒸发和入渗补给精确研究.开始将大气2土壤2地下水作统一体,使不同地下水埋深、不同土质、有无作物覆盖的潜水蒸发影响的研究得以系列化、标准化、促进了水 均衡研究及水资源“三水转化”的评价,但大多数地中渗透仪内无专门的水、盐传感设备,如同“黑匣子”,不能获得其中水分含量( z ,t)盐分含量S或溶质浓度C ( z ,t)的时空分布;受有限设备条件的约束,地中渗透仪的处理不能过多,所设计处理只 能限于有限的观测设备.(4)七十至九十年代,在土壤水盐运动的能量转换理论和地下77 本文于1995年5月30日收到.系世界银行中国北方灌溉资助项目. 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.水数值模拟技术的推动下,在传统地中渗透仪或大型室内土柱内部增加水分或盐分传感 器,进行近代技术的更新改造,使数学模型的数值解与实验技术得以结合,“黑匣子” 得到识解.进而用土壤水盐数据采集系统,对ASGC系统的溶质迁移进行动态模拟.可对复杂环境下进行处理多、效率高、周期短、投资省的精确计算和分析评估;本文属 于此例,以中国河套灌区为目标,研究地下水与土壤盐渍化关系的动态模拟.二、实验方法与测定设备(一)实验物理模型 土柱4根,高220cm、直径14cm,内装均质扰动土(砂土、重砂壤土、轻砂壤土、粉质壤土) .顶部用红外线灯照射模拟稳定蒸发,用标准蒸发皿d= 20cm率定水面蒸发量;底部设排水滤层及地下水位控制装置,模拟不同地下水埋 深;土柱上方设马氏供水箱,模拟降水(灌溉)入渗.构成蒸发2入渗再分配,不同地 下水条件,均质土柱溶质迁移的物理模型.(二)土壤水分数据采集子系统 用有较强垂直空间分辨能力,受土壤化学成分影响较小、易于防护的透射法,放射源为137Cs.射线探头由行走系统操作移动,对土柱轴线不同空间点进行自动采集.(三)土壤盐分数据采集子系统 土柱内固定采集点(间距10cm)预埋接触式盐分传感器,测定土壤溶液电导值EC ,经16通道选通装置自动采集点溶液EC (25)或C值.(四)微机管理系统 经两块数据转换卡将电导仪和射线定标器与PC386计算机连接,用汇编语言编制软件,实现自动选择量程、数据自动采集和记录、存贮的微机 管理. 设土壤含水率目标值为,土壤溶质总浓度目标值为C ,由数据采集系统定时( t)、定点( z)采样获得其时空分布:土壤湿度场记为( Zi,tj) ,溶质浓度场记为C ( Zi,tj) : i= 1,2,I ,I 空间观测点数; j= 1,2, J , J 时间 观测次数;IJ=系统内水分与盐分(浓度)采样品总数.三、数学模型与数值计算方法(一)数学模型 在土壤水2地下水动力学理论指导下,多孔介质中溶质迁移的数学模型据总结7,12,主要是在对流2弥散方程的基础上分别考虑吸附作用,或不动水体的 单独影响.考虑到河套灌区的次生盐渍化特征及河套平原沉积与洪积交错成土特征,土 壤中的盲孔现象(含不动水体)和吸附作用两者均不可忽视,经过专门实验和计算论 证6,10,在对流2弥散方程的基础上需增加同时考虑含有不动水体与吸附作用两种影响的改进模型. 对于均质各向同性土壤,在忽略温度势及生物吸收作用,只研究土壤吸附总量,不 涉及单一离子组分和多离子间相互作用;对于大面积灌排系统中的任一单元体,忽略水 平方向水力梯度影响等假定下,构成一维饱和2非饱和流场,地下水与土壤盐溶质迁移 的数学模型,可分别由水分运动,盐溶质迁移两方程耦合构成:1.土壤水分运动方程 为使数值计算迭代方便,并使计算误差较小和适用于饱和87 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved. 非饱和流,采用以v为因变量的土壤水分运动方程,用隐差分法求解式(1) .A ( v)9v 9t=92v 9Z2+ B ( v)9v 9Z,(11)初始条件: v ( Z , t) = v0( Z) ,Z 0, t =0,(12)上边界条件:一类v ( Z , t) = v0( Z) ,Z =0, t 0,(13)二类-v9v 9Z+ K( v) = R ( t) , Z =0, t 0,(14)蒸发时R ( t) = -Es( t) ,(15)再分配时R ( t) =0,(16)下边界条件9v 9Z=K( v)V,(17)土壤 不式中:Cr为入渗水浓度;Dsh( V ,)=mDsh;Ic( Z ,t)为称源汇项,表示单位时间内、单位体积溶液中由于土壤吸附(解吸)、衰变及沉淀(解析)作用引起的溶质质 量变化.Ic( Z , t) =9S 9t+9Q 9t,(3)式中: S为土壤吸附量.由实验求得含有两种容量因子土壤吸附及解吸过程的河套土壤动力吸附模型9为:S= Smax(1-be-2c) ,(4)式中: Smax为土壤最大吸附量,对轻砂壤土Smax= 9110mmol/ kg;C为溶质平衡浓度;a为经验系数,吸附时a= 010454,解析时a= 010355;b为起始吸附容量因子,吸附 时b= 019505,解析时b= 019765;Q为单位体积土壤中固态盐分含量,采用系数法的下式计算:9Q 9t=j+12i( Ci-Cmax) ,(5)=1, Cj+1 2iCmax, (当j+12i=0时,令=1)2 Cj+1 2i Cmax, 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.流动区与非流动区溶质通过分子扩散交换:9imCim 9t=( Cm-Cim) ,(6)imi= (1-f )i,(7)式中:为质量交换系数;f=m/,f= 0197610. 上述式(2) ,采用单步反向追踪的特征 有限元方法(CFE)求解.联立解式(1) ,(2)即可定解土壤溶质的水盐运动目标.3.模型的正确性 用Naylor和Finger (1977)提出的三步法,在灌溉和冲洗淋 盐、蒸发2入渗再分配条件下进行验证后研究提供了如下简化条件:在蒸发2入渗再分配 时,不动水体具有“汇”的效用,对溶质迁移有较大影响,在长期预测时可忽略吸附作 用的影响.而在降雨或灌溉入渗条件下,土壤吸附作用具有“源”的效应,对溶质迁移有较大影响,因而不动水体影响可忽略.上述成果较Nielsen(1986)总结的单独考虑 吸附作用,单独考虑不动水体影响的数学模型的应用条件作了具体改进补充.(二)数值计算方法 求解土壤水运动v方程比较了隐式差分和中心差分格式,认为虽然中心差分精度较高,但要求步长 t太小,在长期蒸发模拟计算中难以采用,故 选用隐式差分法.求解具有双曲型特征溶质迁移方程,由于浓度峰面太陡,在Peclet数( Pe) ,Courant数( Cr)较大时,用一般数值方法易产生数值弥散(失真)或数值振荡(过 量) ,常在浓度峰面附近出现数值摆动或峰面模糊不清.经对国内外公认的三种先进数 值方法: Bresler时间导数取高次近似的差分;二阶迎风隐式差分;特征2有限元法;经 稳定流、非稳定流的多种典型状态的比较计算认为;由于CFE法将双曲型方程转化为抛物型方程后,消除了数值计算困难的隐患,有精度和效率高的优点,有广泛适应性 如:Pe100,Cr015,故溶质迁移用CFE法6 ,7.四、地下水与土壤盐渍化关系的动态模拟(一)模拟系列设计 根据河套灌区概念模型和灌排条件4,用以下特值划分模拟系列和挡次:1.气候 以水面蒸发强度代表共4个挡次:E0= 110,610(多年平均) ,1010(近似46月平均) ,1815mm/ d.2.土壤 选定灌区内3种代表土类,重砂壤土、轻砂壤土、粉质壤土.土体构型均质扰动土.盐分垂向分布2种,全剖面均匀分布,表层20cm呈现三解型分布,以下 均匀分布.盐分化学组成,以Cl-型为主.积盐动态计算深度5种,Hs= 20,50,100,200,250cm.盐溶质浓度垂向空间分布11种C ( Z)= 0,20,40,60,80,100,120,140,160,180,200cm.3.地下水 埋深6种,Hg= 015,110,115,210,215,310m.矿化度5种,以 单盐NaCl进行实验模拟,Cg= 011,015,110,510,1010g/ L.4.预测时间 3种动态过程,土壤总蒸发时间( tmax) :t1= 1050d、用于两次灌 水间,t2 45d,用于播种 收割,t3= 120d相当全生育期.5.灌溉或冲洗定额 3种,M= 60,80,100m3/亩,灌溉水质Ci= 015110g/ L.08 1995-2005 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.由此组成复杂ASGC系统的多处理、多组合的模拟系统,可表征3种主要自然因 子A、S .G及人工因子M与t的相互关系,即一定条件下地下水与土壤盐渍化动 态7.(二)动态模拟图谱 选择稳定蒸发条