考虑基岩裂隙作用的流域水文过程模拟

石 清，陈 喜，魏玲娜，张志才

(1.河海大学水文水资源与水利工程科学国家重点实验室，江苏南京 210098;2.中国水电顾问集团北京勘测设计研究院，北京 100024;3.南京信息工程大学水文气象学院，江苏 南京 210044)

美国华盛顿大学开发的分布式水文-土壤-植被模型(DHSVM)包含了地表与土壤水流的描述，可以较好地模拟洪水过程［1］。然而，已有研究发现在进行长期水文过程模拟时，现有模型对枯期径流拟合普遍较差［2］。国内新安江模型用于流量过程模拟时，退水段的流量过程主要受地下水消退系数控制，该系数对退水段影响较大［3-4］。在参数调试过程中，退水段的拟合精度有时与整个流量过程的拟合精度相矛盾，故在模拟过程中需综合考虑两者之间的关系，尽量做到兼顾两者的拟合精度［5］。

在我国南方广大山丘区，枯期径流不仅来源于土壤层，还来源于下覆基岩裂隙水。基岩裂隙水是我国分布最广泛的地下水类型之一，浅部裂隙地下水都与上覆土壤水有不同程度的循环交替［6］，进而影响产汇流过程。笔者以东江流域内星丰流域为研究区，基于DHSVM 2.0版本，根据流域地形、土壤、植被覆盖率动态变化特征以及基岩裂隙水对径流的调节作用，建立考虑植被-土壤-基岩裂隙作用的分布式水文模型，模拟红壤丘陵区流域水文过程。

1 流域概况及观测资料

星丰流域地处东江上游，为高丘陵山区，面积42.6 km2，海拔150～508 m，平均比降为3.14%。流域内土壤主要为麻赤红壤、页赤红壤和潴育水稻土，其中麻赤红壤占72%，土层较厚。按国际制土壤分类标准，流域土壤概化为壤土、砂壤土和黏土3种类型。经室内实验分析，土壤密度为1.32～1.45 g/cm3，孔隙率为41% ～45%，土壤含水率为12.3% ～32.7%。

星丰流域位于地下水二级功能区东江河源和平地下水水源涵养区，基岩裂隙水是枯期径流的主要来源［7］。研究区内基岩多属变质岩类裂隙含水岩组，石英砂岩占77.4%，粉砂岩占22.6%。花岗岩风化后形成的风化壳比较深厚，一般在10 m以上，深厚者可达30～40 m。风化壳含有很多石英砂，其黏结力差、透水性强、结构松散［8］。

星丰流域林地占流域面积的91%，耕地占8%，其余为建筑用地和水面。野外实地考察发现，星丰流域现状覆被以松、杉混杂的针叶林为主。基于2001年1月8日至2008年12月18日共184期250 m×250m空间分辨率的MODIS植被指数EVI数据(MOD13Q1)，利用Mu等［9］2007年提出的算法，反演不同类型的植被覆盖率F(图1(a))。基于MODIS数据(MOD13Q1，MOD15A2，MOD09A1)，求得不同类型的植被覆盖率F(图1(a))、典型覆被叶面积指数(L)的月变化曲线(图1(b))以及各种植被类型的短波冠层反射率A(图1(c))。

图1 不同植被的F，L和A月变化过程Fig.1 Monthly variations of fractional coverage F，leaf area index L，and canopy reflection A for different types of vegetation

采用流域内4个雨量站和1个水文站2004—2008年实测逐日降雨和流量观测资料，并移用研究区附近河源气象站气温、风速、相对湿度、长短波辐射等气象资料进行模型计算。

2 DHSVM及改进

2.1 DHSVM结构

DHSVM以地面高程模型(DEM)栅格节点为中心，将流域划分为若干计算栅格单元。根据栅格单元地形特征，描述短波辐射、降雨、气温、坡面流及栅格单元土壤和植被特性时空变化，在每一个计算时段内，根据各栅格的能量和质量平衡方程联立解、栅格单元间坡面流和壤中流汇流演算方法，建立栅格单元水文联系。根据河道分级演算，计算河网汇流过程。DHSVM实现了栅格尺度上流域水文过程及其与下垫面植被、土壤等影响因素的系统描述。

DHSVM中降雨分别被上、下冠层植被截留，直至最大截留能力Icj:

式中:rj——冠层截留系数;Fj——冠层覆盖率。下标j若为o则表示上冠层，若为u则表示下冠层。

冠层截留量变化为

式中:St——t时刻冠层截流量;P——降雨量;EIj——植被潮湿部分蒸发量(即截留雨量蒸发)。

采用双层Penman-Monteith方法模拟上、下冠层散发量，结合土壤水分条件计算土壤实际蒸发量［10］，三部分之和作为实际总蒸散发量。

a.植被潮湿部分蒸发量:

其中

式中:Epj——潜在蒸散发量;Awj——潮湿部分面积，干燥部分则为(1-Awj);Δtw——以Epj蒸发截留雨量所需时间。

b.植被干燥部分散发量:

式中:Ej——植被表面蒸腾量;ζ——饱和水汽压～气温关系曲线的斜率;γ——空气湿度常数;rcj——冠层阻力;raj——空气动力学阻力;Δt——计算时段长。

c.土壤实际蒸发量:

其中

式中Fe为土壤解吸率，由土壤对表层的输水速率确定，是土壤类型和上层土壤含水量的函数。d.土壤层水量平衡采用3层根带模式模拟水分在上、中、下非饱和土壤层中的运动:

式中:dk——第k层土壤的厚度;n——土壤层数;θk——第k层土壤的含水率;If——计算时段内渗入土壤层中的水量;Qv——流入下层的水量;fo，k，fu，k——第k层土壤中对上冠层、下冠层的根系比;Vexk——地下水位上升时地下水对土壤水的补给;Qsin，t，QS，t——计算时段初侧向壤中流的入流量、出流量。

假定土壤饱和渗透系数随土壤深度呈指数递减，推导出土壤导水系数Tx，y的计算公式如下:

式中:Kx，y——栅格(x，y)表层土壤侧向饱和渗透系数;fx，y——指数衰减系数;zx，y——地下水位埋深;Dx，y——土壤层厚度。

采用运动波方法计算侧向饱和壤中流，用Wigmosta［11］发展的准三维路径模式计算坡面流。坡面地表流、壤中流逐栅格汇流并流向流域低洼处汇集进入河道，利用线性槽蓄法进行河道汇流演算。

2.2 模型改进

分析星丰流域2003—2008年实测日流量过程，依据饱和水运动Depuit-Boussinesq方程对退水过程较长的典型年份(2003年、2004年、2006年)退水过程进行拟合(图2)。由图2可知，退水特征明显表现为2个阶段，原因在于南方红壤丘陵区多为二元水文地质结构［12］，退水段径流主要来源于土壤层及基岩裂隙水。上层土壤层渗透性按式(9)进行衰减，退水较快;下层基岩裂隙退水过程较稳定。因此，笔者在饱和土壤层之下添加裂隙含水层模块，综合考虑土壤水和基岩裂隙水之间的水力联系，对DHSVM进行改进。

采用Wigmosta等［13］提出的地形驱动影响下饱和壤中流的计算方法，利用地下水位(而不是地面高程)计算栅格之间的水力梯度。由于水头的变化，需要对每一个计算时段内的地下水位和相应的坡度、坡向、流向、水力梯度进行计算。裂隙含水层的水分运动分为垂向和侧向两部分进行计算，地下水流量的计算公式

图2 星丰流域枯期径流退水过程规律Fig.2 Low-flow curve fitted by Depuit-Boussinesq model for Xingfeng Basin

如下:

式中:fi——栅格(x，y)在i方向的流线宽度;GGW——水力梯度;Kl——裂隙含水层侧向水力传导率;S——栅格单元内含水层含水量;Kd——含水层下基岩的垂向传导率。

裂隙含水层的水量平衡方程为

式中:St，St+Δt——t，t+ Δt时刻裂隙含水层含水量;(QGin，t-QG，t)——裂隙含水层中基流的净增水量(单元栅格中侧向流入量QGin，t减去流出量QG，t);Qv，GW——上覆土壤层向裂隙含水层渗漏的水量;Vex，GW——裂隙含水层水位上升对其上覆土壤层的补给量。

紧邻裂隙层的上覆土壤层水量平衡方程式改写为

当裂隙含水层没有到达饱和时，上覆饱和壤中流以下渗的形式对裂隙含水层进行补给，直至裂隙含水层饱和。当栅格单元内裂隙含水层的含水量大于其最大容量时，多余的水量将“返还”给上覆土壤层。除了和上覆土壤进行水分交换，裂隙含水层水分还会向下渗透，逸出含水层，形成“深层损失”，“深层损失率”由基岩的水力传导率决定。当栅格单元内的地下含水层水位高于河床，裂隙含水层水量将会被河道截留。

3 模型应用

3.1 模型率定与验证

根据DEM将流域划分为50 m×50 m栅格，共162行、156列。模型计算参数率定期为2003—2006年，验证期为2007—2008年，计算步长为2 h。根据模型计算与实测流域出口断面流量过程，采用人工试错法对模型参数进行率定。通过确定性系数和径流深相对误差来评测模型日径流模拟的精度。表1、表2列出了模型主要土壤、基岩裂隙以及植被特征参数的取值。

表1 土壤及基岩裂隙特征参数Table 1 Soil property parameters

表3为利用DHSVM及本文改进的DHSVM进行参数率定与模型验证的结果。如不考虑基岩裂隙水，率定期枯水年份(2003年、2004年)计算径流量R显著偏小，丰水年份(2005年、2006年)则偏大。这是因为DHSVM本身在枯水处理方面有一定缺陷，存在退水衰减过快的问题，枯季径流模拟效果不佳。在模型率定中，由于枯水年径流较小，模拟径流深偏小，为了在连续模拟中保持水量平衡，使丰水年径流模拟结果相应地偏高。进一步分析发现，裂隙含水层对年内径流分配具有一定的调节作用，改进后DHSVM模拟年径流确定性系数有显著提高(表3)，流量过程线的吻合度较高(图3)。

表2 植被特征参数Table 2 Vegetation property parameters

表3 DHSVM改进前后蒸散发及径流模拟结果Table 3 Comparison of evapotranspiration and streamflow before and after DHSVM revision

图3 实测与模拟流量过程线Fig.3 Observed and simulated streamflow

利用实测流量和模型改进前、后的模拟流量序列，采用期望值公式，将流量按从大到小进行排频(图4)，显然，考虑裂隙作用对枯期流量模拟结果改进显著。进一步选取频率为75%的流量值Q75%作为枯季径流的特征值，统计出不同历时下实测和模型改进前、后模拟的枯季流量低于Q75%的发生次数(图5)。通过比较可以发现，改进后枯期径流模拟精度比改进前有较大提高。

3.2 参数灵敏度分析

选用扰动分析法，计算参数在小范围内变化所导致的模型输出变化率，通过敏感度因子Sen加以衡量，计算公式如下:

式中:y，y'——参数变化前、后输出值;x——变化前参数值;Δ——参数变化值。

选取各参数敏感度因子的算数平均值作为该参数的灵敏度值，模型计算输出值选用上覆土壤层含水量、日平均蒸散发量、日平均流量以及汛期日平均流量、枯季日平均流量。

图4 实测与模拟流量历时曲线Fig.4 Observed and simulated flow duration curves

分析表明:对上覆土壤层含水量最敏感的参数是孔隙度Φ和田间持水量θfc，其他参数敏感度都较小;对日平均蒸散发量较敏感的参数有植被参考高度H，F，θfc、最小气孔阻力Rsmin;对日平均径流量最敏感的参数排序为H，F，θfc，Rsmin，其中汛期日径流量敏感性系数较大的参数主要有H，F，Φ和裂隙含水层孔隙度GP，而枯期日径流量敏感性系数较大的参数主要是Φ，θfc，GP。由此可见，上覆土壤含水量及蒸散发主要受土壤、植被参数影响，基岩裂隙含水层参数对汛期和枯期径流起调节作用，减少汛期流量、增加枯期流量，对枯期流量影响显著，但对全年流量影响不大。

图5 各时段实测和模拟流量低于Q75%的数目Fig.5 Occurrences of observed and simulated flow below Q75%

4 结 语

根据东江星丰流域植被、土壤、基岩特征，在考虑裂隙含水层对径流的调节功能情况下对DHSVM进行改进。流域水文过程模拟及参数灵敏性分析表明，基岩裂隙水对径流过程模拟效果影响显著，模型改进后对径流过程的模拟精度有较大提升，尤其对枯期径流的模拟更接近实际观测结果。基岩裂隙含水层参数对汛期和枯期径流起调节作用，减少汛期流量、增加枯期流量，对枯期流量影响尤为显著。

［1］WIGMOSTA M S，VAIL L W，LETTENMAIER D P.A distributed hydrology-vegetation model for complex terrain［J］.Water Resources Research，1994，30(6):1665-1679.

［2］BATHURST JC.Physically-based distributed modelling of an upland catchment using the systeme hydrologique europeen［J］.Journal of Hydrology，1986，87(1/2):79-102.

［3］ZHAO Renjun.The Xin＇anjiang model applied in China［J］.Journal of Hydrology，1992，135(1/2/3/4):371-381.

［4］赵人俊.流域水文模拟:新安江模型与陕北模型［M］.北京:水利电力出版社，1984.

［5］黄国如，陈永勤.基于新安江水文模型的东江流域枯水径流模拟［J］.华南理工大学学报:自然科学版，2006，34(11):93-98.(HUANG Guoru，CHEN Yongqin.Low flow simulation of the Dongjiang Basin based on Xin ＇anjiang watershed hydrological model［J］.Journal of South China University of Technology:Natural Science Edition，2006，34(11):93-98.(in Chinese))

［6］陈宇，温忠辉，束龙仓.基岩裂隙水研究现状与展望［J］.水电能源科学，2010，28(4):62-65.(CHEN Yu，WEN Zhonghui，SHU Longcang.Status and prospect of research on bedrocKfissure water［J］.Water Resources and Power，2010，28(4):62-65)(in Chinese))

［7］广东省水利厅.广东省地下水功能区划［R］.广州:广东省水利厅，2009.

［8］邓南荣，李定强，王继增，等.广东省东江流域土壤侵蚀空间分布特征研究［J］.中国水土保持，1999(5):21-23.(DENGNanrong，LIDingqiang，WANGJizeng，et al.The characteristic of spatial distribution of soil erosion in Dongjiang Basin of Guangdong Province［J］.Soil and Water Conservation in China，1999(5):21-23.(in Chinese))

［9］MU Q，HEINSCH F A，ZHAO M，et al.Development of a global evapotranspiration algorithm based on MODIS and global meteorology data［J］.Remote Sensing of Environment，2007，11(4):519-536.

［10］EAGLESON P.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.

［11］WIGMOSTA M S，PERKINS W A.Simulating the effects of forest roads on watershed hydrology［J］.Water Science and Application，2001，2:127-143.

［12］刘昆，何捷，焦树林，等.粤东花岗岩流域覆盖型风化壳风化特征比较［J］.生态环境，2008，17(5):1937-1941.(LIU Kun，HE Jie，JIAO Shulin，et al.Comparison on the weathering characteristics of the alluvium-covered weathering crust for the granite catchments in eastern Guangdong province［J］.Ecology and Environment，2008，17(5):1937-1941.(in Chinese))

［13］WIGMOSTA M S，LETTENMAIER D P.A comparison of simplified methods for routing topographically driven subsurface flow［J］.Water Resources Research，1999，35(1):255-264.