地下水作用下土壤水蓄量变化及其对蒸发通量影响的模拟

黄远洋,陈　喜,张志才,郑　健,4

(1. 河海大学水文水资源与水利工程科学国家重点实验室,江苏 南京　210098;2. 河海大学水文水资源学院,江苏 南京　210098; 3. 中国科学院重庆绿色智能技术研究院，重庆 400714;4. 浙江省水利水电勘测设计院, 浙江 杭州　310002)



地下水作用下土壤水蓄量变化及其对蒸发通量影响的模拟

黄远洋1,2,3,陈喜1,2,张志才1,2,郑健1,2,4

(1. 河海大学水文水资源与水利工程科学国家重点实验室,江苏 南京210098;2. 河海大学水文水资源学院,江苏 南京210098; 3. 中国科学院重庆绿色智能技术研究院，重庆 400714;4. 浙江省水利水电勘测设计院, 浙江 杭州310002)

摘要：将基于数值积分的土壤水平衡计算与稳态蒸发通量解析解相耦合,提出一种逼近Richards方程的迭代计算模型,分析不同地下水位下土壤水蓄量变化及其对土壤水蒸发量和潜水蒸发量的影响。采用Hydrus-1D模型对迭代模型进行对比验证,并将模型应用于安徽五道沟实验站地下水埋深分别为0.4 m、0.6 m和1.5 m的3个蒸渗仪的土壤水蓄量变化及蒸发通量的模拟。结果表明,迭代模型具有较高的计算精度,能较好地模拟受地下水位以上毛细管力作用的土壤水蓄量动态过程及蒸发通量,还能反演不同地下水埋深下的土壤渗透系数和孔径分布系数。

关键词：Richards方程;土壤水稳定状态;土壤水蓄量;土壤水蒸发;潜水蒸发

土壤蒸发和潜水蒸发通量估算是水文循环的重要环节,也是水资源评价的重要内容。土壤水分运动以及蓄量变化是蒸发通量计算的核心,理论上可采用Richards运动方程[1]推求,如目前常使用的Hydrus-1D模型[2]、FEFLOW模型[3]和SUTRA模型[4]等。该类模型采用有限差分或有限元方法对Richards方程进行求解,容易因时空分辨率选择不恰当而不收敛,且在流域尺度上计算耗时长。采用特定边界条件或简化土壤水动力方程推求水通量的解析解和半经验模型,因计算简单而得到广泛使用。如Gardner[5]、Eagleson[6]以及Famiglietti等[7]给出的土壤水运动上边界基质势无穷大情形下的稳态潜水蒸发量计算公式;Warrick[8]、Salvucci[9]以及Sadeghi等[10]给出的任意稳态蒸发量情况下土壤剖面基质势及含水量分布;雷志栋等[11]和唐海行等[12]基于Gardner公式提出的考虑气象条件潜水蒸发量的半经验公式。但是,该类模型缺乏对土壤水蓄量动态变化的描述。

基于土壤水蓄量动态变化的水平衡计算方法是水文模型构建的基础。但在目前概念性水文模型中,对潜水面以上毛细管力等外应力如何影响土壤水蓄量以及土壤蒸发和潜水蒸发等水文变量还缺乏描述,如新安江模型中的分层蒸发模型[13],在模型概化中对蒸散发作用的土层深度缺乏定量分析,也缺乏显式描述地下水位变动对蒸发通量的影响。基于土壤水运动规律与水平衡计算,Soylu等[14]将根系层的土壤水蓄量作为土壤水流运动上边界条件,改进潜水蒸发量计算的Gardner-Eagleson公式[5-6],并采用数值积分方法对公式进行求解,有效改善了模型的收敛性。但该方法对没有根系层的裸土及根系以下土壤水蓄量推求不适用。

笔者以Sadeghi等[10]提出的稳态蒸发通量解析解及其对土壤水分剖面分布的数学分析为基础,将基于数值积分的土壤水平衡计算与稳态蒸发通量解析解相耦合,提出一种逼近Richards方程的水通量迭代计算模型,使得模型适用于裸土及根系以下土壤,同时具有良好的收敛性。通过对不同地下水埋深下的土壤蒸发和潜水蒸发量的模拟,将迭代计算模型与Hydrus-1D模型进行对比验证。将模型应用于五道沟水文实验站蒸渗仪土壤水蓄量变化及蒸发通量的模拟,反演了不同地下水埋深下的土壤渗透系数和孔径分布系数。

1模型基本原理及计算方法

1.1　土壤水分运动方程

在地下水位以上一维土柱内,非饱和土壤水分的运动可以采用Richards方程进行描述:

(1)

式中:θ——土壤含水率;t——时间;z——以地下水位为参考的垂向高度,向上为正;K——土壤非饱和渗透系数；ψ——土壤基质势。

求解式(1)需要θ和K与ψ的关系,本文采用Brooks-Corey[15]模型:

(2)

其中

Kr=K/Ksp=2+3λSe=Θ/ΘeΘ=θ-θrΘe=θs-θr

式中:Ks——饱和渗透系数;ψb——进气势;λ——土壤孔径分布系数;Se——饱和度;Θe——有效含水率;θs——饱和含水率;θr——残余含水率。

1.2　稳态下地下水面以上土壤水分分布及蓄量表述

稳态下,式(1)中∂θ/∂t=0,即土柱内水通量在垂向上不发生变化:

(3)

地下水位附近的水通量,即稳态潜水蒸发量可写为

(4)

此时,土柱内土壤基质势剖面ψ(z,r)(图1(a)),图中hm为最大干化层厚度)可通过式(2)(4)推求。

令r=Eg/Ks,当r=0时由式(4)可得

(5)

当r>0时,ψ(z,r)可通过其反函数z(r,ψ)进行推求。Sadeghi等[10]推导出式(2)(4)的z(r,ψ)解析式,即当r≤1时:

(6)

其中

当r>1时:

(7)

其中

式中:ψlim——使土壤水分保持连续的土壤基质势极限值。

蒸发作用使土壤基质势增大,土壤含水率减小。当近地表土壤基质势超过ψlim后,裸土表面发生干化现象[5, 16-18]。干化层内的ψlim接近大气势 (图1)：

图1　稳定状态下的土壤剖面Fig. 1　Soil profile in steady state

其中

h=D-zlim

(8)

式中:ψatm——大气势,1 000 m时可代表一般大气条件;D——地下水埋深;h——干化层厚度;zlim——毛管作用的最大高度。

Sadeghi等[10]按干化层底部土壤基质势趋于无穷大推导出zlim:

(9)

将基于式(5)(6)和(8)推求的土壤基质势剖面ψ(z,r)代入式(2),即可推求对应于某一r时的土壤水分剖面Θ(z,r)(图1(b))。将Θ(z,r)在土柱内进行数值积分即可推求土壤水蓄量W(r):

(10)

1.3　非稳态下土壤水蓄量变化及蒸发计算

1.3.1潜水蒸发量与土壤水蓄量的耦合方法

一维土柱内的水分运动处于非稳定状态时,土柱在时段长度△t内的水量平衡方程由式(1)可得

(11)

式中:W1、W2——土柱在时段初和时段末的土壤水蓄量。

无雨期式(11)右边第一项为上边界z=D处的通量(即△t内的土壤蒸发量Es),右边第二项为下边界z=0处的通量(即△t内的潜水蒸发量Eg),则式(11)可写为

(12)

Es可采用一层蒸发模式[13, 19]计算:

(13)

其中

Ep=ηE601

式中:Ep——潜在蒸发量;E601——E601蒸发皿蒸发量;η——蒸发皿-潜在蒸发转换系数;Wmax、Wmin——最大、最小土壤水蓄量。

由式(12)可知,土壤水蓄量的变化受Es的消耗作用和Eg的补充作用影响。如在毛细管力作用下,Eg接近Es,则土柱内土壤水蓄量近似稳定状态[14, 16]。因此,采用稳态土壤基质势剖面(式(5)～(9)),由式(10)推求W2:

(14)

在极端干旱气候条件以及特定的土壤特性和地下水位条件下,土水势梯度最大,干化层厚度达最大hm(由Eg,max通过h=D-zlim推求),此时土壤水蓄量达到最小:

(15)

式中:Eg,max——相应于hm的最大潜水蒸发量,通常可以根据观测数据推求。

当蒸发通量为零时(即Eg=0,r=0),土水势梯度最小,土壤水蓄量达到最大:

(16)

将式(13)～(16)及由迭代方法求解的Eg代入式(12)可推求土壤水蓄量的变化。

1.3.2迭代求解方法

2模型验证与应用

2.1　模型验证

构建如下一维土柱蒸发算例:上边界条件为土壤蒸发速率,假设其逐日过程符合均值为1.7 mm、方差为2 mm2的Gamma分布[20];土柱下边界条件为地下水位,埋深分别为0.4 m、0.6 m和1.5 m。考虑3种土壤质地:砂质壤土、砂质黏壤土和粉质黏壤土,基于Rawls等[21]给定的3种土壤的Brooks-Corey模型参数,分别采用Hydrus-1D模型和本文基于蓄量平衡的水通量计算迭代模型,对上述9种情形下的潜水蒸发量(即下边界水通量)及土壤水蓄量变化进行模拟,结果如图2所示。

由图2可知,9种情形下Hydrus-1D模型和迭代模型所计算的土柱日潜水蒸发量(图2(a)~(c))和日蓄量变化量(图2(d)～(f))的散点都均匀分布在1∶1线附近;当地下水埋深较大时,迭代模型的结果相对偏大,这是由于稳定状态假定在地下水埋深较大相对较难达到。总体而言,迭代模型具有较高的计算精度。

2.2　模型应用

将基于蓄量平衡的迭代模型应用于安徽五道沟大型水文水资源综合实验站3个蒸渗仪的潜水蒸发过程的模拟和分析,地下水埋深分别为0.4 m、0.6 m和1.5 m,装填原状砂姜黑土,表层无作物,面积均为0.3 m2。数据选取1992—1994年和1996—1999年降雨量、E601蒸发皿蒸发量、3个蒸渗仪的潜水蒸发量、入渗补给量和地表径流量等观测资料。将降雨结束的第三日作为无雨期的开始,假定此时的蒸渗仪水分剖面近似于通量为零时的稳态剖面,即式(5)和式(10)所描述的剖面,计算土壤含水量及潜水蒸发量。

根据五道沟实验站土壤物理性质分析资料[22],土壤有效含水率Θe=0.3。通过对E601蒸发皿蒸发量和地下水埋深为0 m的蒸渗仪潜水蒸发量进行回归分析,给定η=1。

图2　Hydrus-1D模型和迭代算法计算的日潜水蒸发和日蓄量变化量散点图Fig. 2　Scatter plots of daily groundwater evaporation and water storage simulated by Hydrus-1D and iterative method

2.2.1模型参数率定与验证

分别采用1992—1994年和1996—1999年的潜水蒸发量进行参数率定和模型验证。以Eg的相对误差δ作为目标函数,采用SCE-UA[23]优化方法分别率定3个不同地下水埋深蒸渗仪内的土壤参数。初步率定后给定ψb=10 cm,然后对Ks和λ进行优化,参数率定结果和δ见表1。

表1　率定的模型参数和目标函数值Table 1　Calibrated values of model parameter and objective function

由表1可知,Ks符合在五道沟实验站内对原状砂姜黑土采用单环入渗法[24]分层测定的结果,即Ks在10～80 cm深度内(不计0～10 cm的耕作层)为9.6～163 mm/d;λ取值在Rawls等[21]给出的砂质黏壤土的λ范围内。

从潜水蒸发量模拟效果来看(表1),模型计算与实测总量的相对误差δ均小于10%。从潜水蒸发系数αEg(即Eg/Ep)来看(表2),模型模拟值与实测值在率定期与验证期都十分接近。

表2　模型计算与实测的潜水蒸发系数Table 2　Simulated and observed groundwater evaporation coefficients

2.2.2地下水位变化对土壤缺水量及蒸发量的影响

地下水位变动影响饱和-非饱和界面水通量如潜水蒸发量,进而影响土壤蓄水量和土壤蒸发量。模拟结果显示(表3),在模拟期和验证期平均潜在蒸发率2.54 mm/d情形下,当地下水埋深从0.4 m增加到1.5 m时,干化层厚度hm从15 cm增到大20 cm,平均潜水蒸发量Eg从0.63 mm/d减小到0.02 mm/d,平均潜水蒸发系数αEg从0.248减小到0.008,土壤蒸发量Es从1.73 mm/d减小到1.52 mm/d,最大蓄水容量(WM=Wmax-Wmin)从21.0 mm增大到46.4 mm。上述结果表明:随着地下水埋深的增大,Eg和αEg显著减少,蓄水量减小,Es减小,缺水量增加,hm增大,且在埋深从0.4 m增加到0.6 m时,土壤和潜水蒸发通量变化显著。

表3　不同地下水埋深下模型计算的土壤水蒸发量、潜水蒸发量及潜水蒸发系数Table 3　Soil water evaporation, groundwater evaporation, and groundwater evaporation coefficient at different groundwater depths simulated using iterative model

3结论与展望

a. 将基于数值积分的土壤水平衡计算方法与稳态蒸发通量解析解相耦合的迭代模型,具有较高的计算精度, 能较为可靠地描述地下水位对蒸发过程的影响机理。

b. 地下水位埋深0.4 m、0.6 m和1.5 m下,砂姜黑土的饱和渗透系数分别为61 mm/d(7.1×10-7m/s)、20 mm/d(2.3×10-7m/s)和34 mm/d(3.9×10-7m/s),孔径分布参数分别为0.58、0.49和0.29,它们分别反映了地下水面以上土壤的渗透性和黏粒含量的平均特征。

c. 地下水埋深0.4 m、0.6 m和1.5 m下,由蒸发作用引起的砂姜黑土潜水蒸发量分别为0.63 mm、0.08 mm和0.02 mm,干化层厚度分别为15 cm、20 cm和20 cm,土壤缺水量分别为21 mm、29 mm和47 mm。

本文提出的迭代模型,适用于地下水作用显著的裸土和根系层以下土壤,参数物理意义明确,数值方法收敛稳定。但对层状土壤的适用性还有待进一步讨论。

参考文献：

[1] RICHARDS L A. Capillary conduction of liquids through porous mediums[J]. Physics, 1931,1(5):318-333.

[3] TREFRY M G, MUFFELS C. FEFLOW: a finite-element ground water flow and transport modeling tool[J]. Ground Water, 2007,45(5):525-528.

[4] VOSS C I, PROVOST A M. SUTRA, a model for saturated-unsaturated variable-density ground-water flow with solute or energy transport[R]. Reston: US Geological Survey, 2002.

[5] GARDNER W R. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[J]. Soil Science, 1958,85(4):228-232.

[6] EAGLESON P S. Climate, soil, and vegetation 3: a simplified model of soil moisture movement in the liquid phase[J]. Water Resources Research, 1978,14(5):722-730.

[7] FAMIGLIETTI J S, WOOD E F. Multiscale modeling of spatially variable water and energy balance processes[J]. Water Resources Research, 1994,30(11):3061-3078.

[8] WARRICK A W. Additional solutions for steady-state evaporation from a shallow water table[J]. Soil Science, 1988,146(2):63-66.

[9] SALVUCCI G D. An approximate solution for steady vertical flux of moisture through an unsaturated homogeneous soil[J]. Water Resources Research, 1993,29(11):3749-3753.

[10] SADEGHI M, SHOKRI N, JONES S B. A novel analytical solution to steady-state evaporation from porous media[J]. Water Resources Research, 2012,48(9):W9516.

[11] 雷志栋, 杨诗秀, 谢森传. 潜水稳定蒸发的分析与经验公式[J]. 水利学报, 1984,15(8):60-64. (LEI Zhidong, YANG Shixiu, XIE Senchuan. The analysis and the empirical formula for the steady-state groundwater evaporation [J]. Journal of Hydraulic Engineering, 1984,15(8):60-64. (in Chinese))

[12] 唐海行, 苏逸深, 张和平. 潜水蒸发的实验研究及其经验公式的改进[J]. 水利学报, 1989,20(10):37-44. (TANG haixing, SU Yishen, ZHANG Heping.The experimental study and the improvement of empirical formula for the groundwater evaporation [J]. Journal of Hydraulic Engineering, 1989,20(10):37-44. (in Chinese))

[13] 赵人俊. 流域水文模拟:新安江模型与陕北模型[M]. 北京: 水利电力出版社, 1984.

[14] SOYLU M E, ISTANBULLUOGLU E, LENTERS J D, et al. Quantifying the impact of groundwater depth on evapotranspiration in a semi-arid grassland region[J]. Hydrology and Earth System Scineces, 2011,15(3):787-806.

[15] BROOKS R H, COREY A T. Hydraulic properties of porous media[M]. Fort Collins: Colorado State University, 1964.

[16] GOWING J W, KONUKCU F, ROSE D A. Evaporative flux from a shallow watertable: the influence of a vapour-liquid phase transition[J]. Journal of Hydrology, 2006,321(1/2/3/4):77-89.

[17] BITTELLI M, VENTURA F, CAMPBELL G S, et al. Coupling of heat, water vapor, and liquid water fluxes to compute evaporation in bare soils[J]. Journal of Hydrology, 2008,362(3/4):191-205.

[18] GRAN M, CARRERA J, MASSANA J, et al. Dynamics of water vapor flux and water separation processes during evaporation from a salty dry soil[J]. Journal of Hydrology, 2011,396(3/4):215-220.

[19] ALLEN R G, PEREIRA L S, RAES D, et al. Crop evapotranspiration: Guidelines for computing crop water requirements[M]. Rome: FAO, Water Resources, Development and Management Service, 1998.

[20] RYSZKOWSKI L. Landscape ecology in agroecosystems management[M]. Boca Raton: CRC Press, 2001.

[21] RAWLS W J, BRAKENSIEK D L, SAXTON K E. Estimation of soil water properties[J]. Trans Asae, 1982,25(5):1316-1320.

[22] 王振龙, 章启兵, 李瑞. 淮北平原区水文实验研究[M]. 合肥: 中国科学技术大学出版社, 2011.

[23] DUAN Q, SOROOSHIAN S, GUPTA V. Effective and efficient global optimization for conceptual rainfall-runoff models[J]. Water Resources Research, 1992,28(4):1015-1031.

[24] CHENG Qinbo, CHEN Xi, CHEN Xunhong, et al. Water infiltration underneath single-ring permeameters and hydraulic conductivity determination[J]. Journal of Hydrology, 2011,398(1/2):135-143.

Modeling of variation of soil water storage with groundwater table and

its influences on transpiration flux

HUANG Yuanyang1,2,3, CHEN Xi1,2, ZHANG Zhicai1, 2, ZHENG Jian1, 2, 4

(1.StateKeyLaboratoryofHydrology-WaterResourcesandHydraulicEngineering,

HohaiUniversity,Nanjing210098,China;

2.CollegeofHydrologyandWaterResources,HohaiUniversity,Nanjing210098,China;

3.ChongqingInstituteofGreenandIntelligentTechnology,ChineseAcademyofSciences,

Chongqing400714,China;

4.ZhejiangDesignInstituteofWaterConservancyandHydroelectricPower,Hangzhou310002,China)

Abstract:An iterative model was proposed to approximate the Richards’ equation by coupling the soil water budget method and the analytical solution of steady-state transpiration flux, and was used to analyze the variation of soil water storage with groundwater table and its influences on soil water evaporation and groundwater evaporation. The coupled model was verified against the Hydrus-1D model. Then, the model was applied to simulation of the variation of soil water storage and transpiration flux at the groundwater depths of 0.4 m, 0.6 m, and 1.5 m at the Wudaogou Hydrological Experimental Station in Anhui Province, and the results were compared with observations from three lysimeters. The results indicate that the coupled model can effectively simulate the effects of the groundwater table on soil water storage and evaporation flux under capillary forces, and conduct a reasonable inversion of the hydraulic conductivity and pore size distribution index at different groundwater depths.

Key words:Richards’ equation; steady state of soil water; soil water storage; soil water evaporation; groundwater evaporation

中图分类号：P339

文献标志码：A

文章编号：1000-1980(2015)06-0562-07

通讯作者：陈喜，教授。E-mail: xichen@hhu.edu.cn

作者简介：黄远洋(1985—)，男，重庆人，博士研究生，主要从事水文物理规律模拟研究。E-mail: huangyuanyang@cigit.ac.cn

基金项目：国家自然科学基金(51079038，40930635，51190091)；江苏省普通高校研究生科研创新计划(CXLX12��0250)

收稿日期：2014-09-26