地球陆面地形地貌复杂多样,既有高低起伏的山地、丘陵,又有地势平坦的冲积平原,构成了诸多山丘与平原交错分布的地理特征。水文循环过程以流域为单元,径流从坡面汇至河道,从上游陡峭的山丘区沿水系逐级汇至平缓的平原区,整个过程始终发生着地表水与地下水的相互作用与转化^[1]。但是由于地形地貌的显著差异,这种相互作用与转化关系存在很大的不同,因而地表水与地下水耦合模拟方式也应有所区别。

地下水过程是流域水文循环过程的重要环节,地下水模拟则是流域水文循环模拟不可缺少的重要组成部分^[2],从这个意义上讲,任何水文循环模型都应是地表水地下水耦合模型,区别在于两者的耦合方式和模拟方法。目前被广泛接受的一种耦合方式分类为紧密型、松散型以及介于两者之间的半松散型^[3]。这种分类方式也与地表水、地下水的模拟方法有一定关系。例如,模型用Saint-Venant方程模拟坡面流,用Richards方程模拟地下水流,两者均采用数值解法,可以实现地表水地下水模拟的同步求解与实时交互,即所谓的紧密耦合模型,如CATHY、InHM、HydroGeoSphere、OpenGeoSys、ParFlow、MIKE-SHE、MODHMS等^{[4,5,6,7,8,9]}。这种耦合模型物理机制强,但数据要求高,模型计算效率低,不适合大时空尺度的水循环模拟。模型地表水文采用集总式或半分布式模拟,地下水采用线性水库模拟,两者之间的耦合仅在于通过信息交互实现各自时段末的水量平衡演算,形成一种相对独立又相互关联的松散耦合方式,如SWAT、VIC、TOPMODLE、三水源新安江模型等^{[10,11,12,13,14]}。这种耦合模型物理机制相对较弱,但数据要求相对较低,且没有迭代计算,模型运行效率高,适合大型流域水文循环的长时间系列模拟。为了使模型尽可能兼具强物理机制和高运行效率,充分汲取紧密耦合模型与松散耦合模型的优势,根据地表水文特点和地下水流模拟需求,将地下水数值模型耦合到半分布式水文模型中,形成地表水地下水的半松散耦合方式,如SWATMOD、GSFLOW、TOPMODEL-MODFLOW、SWAT-MODFLOW等^{[15,16,17,18]}。这种耦合方式在充分考虑水文条件复杂性的同时减少了坡面流数值计算的工作量,还能模拟地下水的侧向流动和地下水位动态,是目前地表水地下水耦合模拟研究的热点之一,也是较有实践应用前景的模型耦合方式。

由于大型山丘平原交错区(如面积在10万km²以上区域)地形地貌和水文地质条件的复杂性,紧密耦合模型在这类区域的使用仍然有一定的局限性^[16];平原区是社会经济活动主要聚集区,水资源开发利用强度大,在水资源匮乏地区地下水过度开采,有些地区甚至产生了严重的地下水降落漏斗,对区域生态环境和可持续发展产生了不良影响^[19,20],需要在区域水文循环整体中开展地下水流精细化模拟研究,显然松散耦合模型难以胜任这项工作。半松散耦合模型为大型山丘平原交错区的地表水地下水耦合模拟提供了一种更加可行的方案。为了提高模拟效率,模型甚至可以仅对平原区地下水开展数值模拟,不同地形地貌区采用不同的耦合模拟方式。

本文将介绍一种基于子流域划分的半分布式水文模型^[11]与基于有限差分解法的地下水数值模型的耦合方式,形成的地表水地下水耦合模型可以仅在平原区(或者关键区)开展地下水数值模拟,平原区地表径流以及山丘区地表径流、地下水流等水文过程仍然采用半分布式模拟。模型将在中国东北地区的三江平原(面积为10.57万km²)进行水文循环模拟应用。该地区山丘区与平原区交错分布,地表水、地下水、人工取用水、区域外调水等相互作用,部分区域地下水超采严重,对当地生态环境和粮食生产安全产生了重要影响^[21]。亟需开展基于全区域水文循环模拟的地表水地下水相互作用与转化研究,具体包括：① 三江平原地表水地下水耦合模型构建与校验;② 全区域水文循环转化关系量化分析;③ 平原区地下水动态分析;④ 山丘和平原分区地下水补给排泄关键参数分析。研究成果可为三江平原水文水资源相关研究与实践提供参考。

目前已公开的地表水地下水半松散耦合模型一般并非原生模型,大都是对已有水文模型与地下水数值模型的改良式耦合^{[15,16,17,18]}。其中,以子流域或水文响应单元(HRU)划分的半分布式水文模型与基于网格单元划分的地下水数值模型的耦合最为典型,相关的研究和应用也比较多^{[15-16,18,22-30]}。这些半松散耦合模型在耦合机制上主要面临2个方面的问题：一是建立2个模型的时空对应关系,二是实现2个模型之间的信息交互。不管是两个模型在程序代码上完全整合,形成一个高度集成的综合性水文模型^[30],还是通过额外的程序智能调用2个模型,实现地表水地下水联合模拟^[27],抑或人工手动建立2个模型之间的数据文件通讯,完成相应水文过程的模拟^[23,24],都需要解决上述2个方面的问题。

半分布式水文模型在模拟流域/区域水文过程时,需要将流域/区域进行细化,如划分为若干子流域;为了体现子流域内下垫面的非均一性对水文过程的影响,又按照土地利用、土壤、坡度、土地管理类型等进一步细分为水文响应单元(HRU)^[11,16,31];而地下水数值模型则需要按照地下水流偏微分方程的数值解法将含水层进行离散,如按照有限差分解法将含水层剖分为网格单元^[32]。2个模型在空间离散上存在很大差异,要实现两者的耦合,首先要建立不同模拟单元的空间叠加和从属关系,为信息交互奠定基础。目前比较流行的做法是建立HRU-网格单元(Grid Cell)、主河道(Main Channel)—河流单元(River Cell)之间的空间映射关系^[26,30]。在模拟的时间步长上,半分布式水文模型一般以日为模拟时段,模拟结果可以按日、月、年、多年平均等不同尺度统计输出^[31];而地下水数值模型可以根据需求设置不同的模拟时间步长,由于地下水流相对地表水流运动更加缓慢,时间步长往往设置得比水文模型更大,因此在模型耦合时需要建立两者不同时间步长的变换关系。为便于数据处理和结果统计,大多数耦合模型在使用时一般都随水文模型的时间离散方式以日为模拟步长^[16,27]。耦合模型在已建立的空间从属关系和时间步长对应关系的基础上,按照既定的计算顺序,依次模拟地表水文和地下水流过程,期间2个模型可以实时通讯(逐时段信息交互),也可以定期通讯(需要的时候再传递数据),通讯的内容可以包括降水/灌溉深层渗漏、河流/水库/湖泊/湿地渗漏、地下水开采量、潜水蒸发、基流、地下水位等。

从已有半松散耦合模型的研发和应用来看,大部分采用HRU—网格单元的空间关联和信息转换机制,有的模型甚至将组成HRU的栅格(如位图的像素)作为地下水模型的网格单元^[23,26],这很大程度上限制了模型在大型区域的使用,不能充分发挥半松散耦合模型在大型流域/区域水文循环模拟中的优势。另外,为了提高模拟效率,半分布式水文模型在划分HRU时,往往会将面积比例较小的HRU(非主要HRU)删除并合并到其他HRU上(用户自行设定阈值),这会导致HRU丢失空间位置属性^[33],这与半松散耦合方式相冲突,因为地下水网格单元需要与空间确定性的HRU建立关联。建立在完全HRU基础上的耦合模型会大幅增加模拟负担,也不利于耦合模型在大型流域/区域的应用。

为了使半分布式水文模型与地下水数值模型能够更加灵活、便捷地耦合,减少建模数据处理工作量和模型运行负担,提高模型在大型流域/区域水循环模拟中的适用性,本文介绍一种新的耦合机制,将三维有限差分地下水流模型整合到半分布式水文模型MODCYCLE^[11]中,形成一套综合性的地表水地下水耦合模型,并将其应用于面积10万km²以上的大型山丘平原交错区,突出展示半松散耦合模型的优势。

MODCYCLE模型研发之初主要面向人类活动影响下的“四水转化”模拟,地表水文采用子流域划分的半分布式模拟,地下水过程采用线性水库模拟^[11]。为了提高地下水模拟的物理性和精度,增加地下水位、地下水侧向流动等模拟功能,模型引入了三维有限差分地下水流模型^[32],并在模型代码中将其作为一个地下水模块嵌入其中,形成一套高度集成的综合性水文模型。

MODCYCLE模型采用C++语言编程,以Access数据库为模型的数据输入输出平台,面向对象模块化的程序设计理念提高了模型的可扩展性,全方位的并行运算技术提高了模型的运行效率。在模型整体结构上,首先利用数字高程模型(DEM)将流域/区域(一个区域可能包含多个流域)划分为不同的子流域。每个子流域对应一条主河道,各子流域之间的水力联系由主河道的树状拓扑结构建立。根据土地利用、土壤类型和土地管理模式的空间叠加,将每个子流域进一步划分为若干HRU。HRU不等同于一块土地,而是具有相同土地利用、土壤类型和土地管理模式的集合体,在子流域内不具有空间位置属性,各HRU相对独立,没有水力联系。每个子流域的河流系统由一条子河道和一条主河道组成。子河道收集来自HRU的水量,部分被输送到子流域内的天然水域(池塘、湿地、湖泊等),部分被输送到主河道。主河道和水库共同构成了模型的河网系统。

模型在子流域尺度将地下水系统分为浅层和深层,但不同子流域地下水是相互独立的,不存在横向地下水交换。这种地下水模拟方式在山丘区可能是合理的,因为山丘区地表水流域通常与地下水流域是重叠的(喀斯特山区除外)^[34,35]。但在平原区子流域分水岭不明显,不同子流域地下水含水层相互联系,地下水横向流动不容忽视,尤其在地下水开采剧烈的平原区。为此,MODCYCLE在原有水文循环模拟框架的基础上,将地下水数值模型纳入其中,在实际应用时可以将大块、完整的平原区从流域/区域分离出来,专门开展地下水数值模拟。

在大型流域/区域开展长期水文循环的高效率模拟,同时兼顾重点区域、重点水文过程的精细化模拟,是MODCYCLE模型的重要优势之一,其技术要点在于地表水地下水模拟的耦合机制。该模型的耦合也需要解决前述2个方面问题：建立时空对应关系和实现信息交互,但又不同于其他半松散耦合模型的解决方法。

2.2.1 时空关联

MODCYCLE模型的地表水文部分采用日步长模拟。为了便于数据信息的实时交互,地下水模拟部分也统一采用日步长,但可根据需求输出不同时间分辨率的模拟结果,由此建立了地表水地下水耦合模拟的时间关联。本节重点介绍模型的空间关联机制。

鉴于前述HRU—网格单元空间关联的弱点,MODCYCLE模型摒弃了这种关联方式,而是采用了子流域—网格单元空间关联的新方式,减少了大量HRU相关数据前处理的工作量,也不影响非主要HRU的精简合并,在很大程度上进一步扩大了半松散耦合模型在大型流域/区域水文循环模拟中的优势。模型子流域和网格单元空间叠加和对应关系见示意图1。

通过子流域与网格单元的空间叠加,建立两者的对应关系,即空间关联。图1b为2种空间剖分格式叠加后的效果。在平原区含水层进行网格剖分时,有些网格单元没有对应的子流域,称为无效网格单元;反之,称为有效网格单元。对于某个子流域而言,所有有效的网格单元都能与之建立关联,区别仅在于两者的叠加关系。以网格单元与该子流域的叠加面积和该子流域的面积之比建立两者的空间关联：

式中： Afr(i,j)为编号为j的网格单元占编号为i的子流域的面积比例; A(i,j)为编号为i的子流域与编号为j的网格单元的叠加面积; A(i)为编号为i的子流域的面积。有叠加关系时, Afr(i,j)>0;反之, Afr(i,j)=0。

2.2.2 信息交互

半分布式水文模型与地下水数值模型通过相互传递必要的数据信息实现两者的耦合模拟。其中,地下水数值模型所必需的源汇项数据,如降水/灌溉入渗补给、地表水体渗漏补给、地下水开采、基流、潜水蒸发等,可以从水文模型模拟获得,而且有些源汇项数据比单独的地下水数值模型(如MODFLOW)输入的源汇项数据更具合理性。水文模型考虑了土地利用、土壤类型、土地管理模式的空间差异对降水入渗的影响;考虑了气候气象条件、农作物类型和需水规律、土壤属性等因素对农业灌溉地下水开采和井灌回归的影响;考虑了气候条件、植被类型及空间分布差异对潜水蒸发的影响;还考虑了各源汇项与其他水文环节的相互作用,这些都是单独的地下水数值模型在处理源汇项数据时所难以达到的。然而,因为地下水模拟的集总性,使水文模型在模拟地下水位/埋深方面存在很大限制,而地下水位/埋深又是模拟地表水体渗漏、基流排泄、潜水蒸发等水文过程的必要数据。耦合了地下水数值模型的水文模型较好地弥补了这一缺陷,即通过信息交互将地下水数值模型模拟的地下水位传递给水文模型,辅助其完成相关水文过程的模拟,包括前述地下水模拟所需的各项源汇项。

MODCYCLE模型基于建立的子流域—网格单元空间关联,以日尺度为模拟步长,开展水文模拟过程中的实时信息交互。模拟开始后模型首先利用给定的初始条件模拟第一日的水文循环通量和状态量,当日模拟结束时将有关通量传递给地下水数值模拟模块,这些通量包括补给项和排泄项2类。从地表水文部分的子流域传递给地下水数值模拟部分的网格单元,该过程根据两者的空间关联进行,计算公式如下：

式中：Q_cel(j)为编号为j的网格单元的源汇强度,为补给量与排泄量的差值,正值表示该网格单元为净补给,负值为净排泄;R_sub(i)为编号为i的子流域向含水层补给的通量,包括降水入渗、地表水体渗漏、灌溉渗漏等补给项;D_sub(i)为编号为i的子流域从含水层排泄的通量,包括潜水蒸发量、基流、地下水开采量等排泄项;I为子流域总数。该公式的意义表示,所有与网格单元j有叠加关系的子流域,将其补给和排泄量按照空间关联的比例分割后汇总到网格单元j,形成该网格单元的源汇项总量。

地下水数值模拟模块利用地表水文模拟传递的源汇项强度,结合给定的边界条件、水文地质参数、初始水头等,模拟当日结束时各网格单元的地下水位/埋深,并利用子流域—网格单元空间关联转换为子流域的地下水位/埋深,协助模型次日的地表水文模拟：

式中：H_sub(i)为编号为i的子流域的地下水位/埋深;H_cel(j)为编号为j的网格单元的地下水位/埋深;J为有效网格单元总数。该公式的意义在于,所有与子流域i有叠加关系的网格单元,将其地下水位/埋深按照空间关联的叠加面积加权平均到子流域i,形成该子流域的平均地下水位/埋深。

模型以子流域—网格单元空间关联为基础,逐日转换、传递模拟的补给排泄量和地下水位/埋深,完成每个时间步长(日尺度)的地下水数值模拟和地表水文模拟,并根据应用需求输出各种时间分辨率的水文循环通量和状态量。

为了验证MODCYCLE模型在大型山丘平原交错区的适用性,选取三江平原作为研究区,开展模型的构建工作,主要包括数据处理和模型校验。

三江平原位于黑龙江省东部,黑龙江及其支流松花江、乌苏里江交汇的区域,由三江低平原、兴凯湖平原、山间河谷平原及交错分布的小兴安岭、完达山、长白山余脉等山丘区组成,面积10.57万km²,属于典型的大型山丘平原交错区(图2)。其中,三江低平原和兴凯湖平原是主要的平原区,社会经济活动强烈,湿地退化严重,井灌水稻发展迅速,部分区域地下水超采严重,亟需开展以区域水文循环模拟为基础的地表水地下水相互作用与转化研究,为三江平原水资源合理开发利用、地下水保护、湿地修复、粮食生产安全提供科学依据。

MODCYCLE作为半松散耦合模型,其地表水文采用半分布式模拟,利用DEM和实际河网划分子流域、生成主河道是模型构建的基础。为了精细化模拟平原区地下水的动态,揭示平原区人类活动影响下地表水地下水的相互作用与转化机制,将全区子流域分割为平原区子流域和山丘区子流域,其中对平原区地下水进行数值模拟。对于跨越山丘区和平原区的子流域,利用既定的山丘平原分界线强制分割,对应的主河道也在分界线处断开,保持子流域与主河道的一一对应关系。全区共划分为1705个子流域,单个子流域平均面积约62 km²(图3)。模型仅在三江低平原和兴凯湖平原开展地下水数值模拟,将含水层垂向分为深、浅2层,水平方向以2 km的间距剖分为206行206列,单层计算单元42436个,其中平原区范围内的有效单元13076个(图4)。

模型构建所需数据除了上述基础空间数据,还需要土地利用、土壤、农业种植结构等空间分布数据;水文地质参数、边界条件、初始条件等地下水数值模拟所需数据;降水、气温、风速、辐射、空气湿度等气象驱动数据;农业灌溉、工业/生活/生态用水、退水、水库调蓄、区域外调水、地下水开采等人类驱动数据。地下水数值模拟部分无需输入源汇项参数,而是由地表水文模拟部分提供,但增加了子流域与网格单元的空间关联数据。

根据模型所需数据的收集情况,本次模拟构建了三江平原2000—2014年水文循环模型。这一时期正是三江平原水稻面积急剧扩张的关键期,地下水开采量骤增,显著干扰了当地的自然水文循环过程,具有很强的代表性,因此有开展模拟研究的必要性。

数据库建立后,模型能够运行,但由于初始变量和参数设置不一定符合实际,使模型的模拟效果一般不能满足需求,因此需要重新调整相关参数,提高模型刻画实际水文过程的能力,即所谓的模型参数率定。为了检验模型预测未来水文过程的能力,还需要用实测数据对率定后的模型进行验证。重复上述过程,直到模型达到符合某个标准的效果。

根据实测数据的收集情况,本文采用2000—2008年的数据进行参数率定,用2009—2014年的数据进行模型验证,校验过程采用人工手动调参方式。模型校验效果采用纳什效率系数(NSE)和决定系数(R²)评估：

式中：Q_{m, i}为第i个流量实测值;Q_s,i为第i个流量模拟值; Q̅m为长系列流量实测值的平均值;O_i为第i个实测值; O̅为长系列实测值的平均值;P_i为第i个模拟值; P̅为长系列模拟值的平均值;n为系列长度或实测值数量。

NSE用于评估模拟径流量与实测径流量的拟合效果,R²用于评估径流量、地下水埋深模拟值与实测值的拟合效果。

本次模型校验共收集了研究区内6个代表性水文站和平原区21个代表性观测井的流量、地下水位长系列实测数据,分别用于校验模型在模拟地表径流和地下水动态方面的效果。图5显示了部分模拟值与实测值的拟合效果,其中,图5a和5b分别为乌苏里江重要支流挠力河宝清水文站(S1)和松花江支流梧桐河宝泉岭水文站(S2)月径流量过程的拟合效果,根据Moriasi的评估标准^[36],NSE和R²均达到了满意标准以上。图6a和6b为三江低平原观测井地下水埋深变化过程的拟合效果,模拟水位变化周期和幅度基本与实测水位过程相符合。图6c为兴凯湖平原观测井拟合效果,受多种因素影响,地下水埋深周期性规律不明显,拟合效果相对较差,但水位变化幅度较小,且模拟值与实测值变化处于同一水平,基本满足地下水模拟需求。图5c和6d分别为所有水文站和观测井全模拟期的模拟值与实测值的对比效果,总体上模拟效果良好。表1列出了全部代表性观测井和水文站模型校验后的NSE和R²。

模型校验过程中,模拟结果对有些参数较为敏感,参数较小的调整能引起地表径流、地下水位较大的变化,这是导致水文循环模拟不确定性的重要因素之一^[37]。根据以往建模的经验,通过反复试算,选取径流和地下水响应敏感的几个关键参数,并不断缩小这些参数的取值范围,以确定模型的最佳参数组合。三江平原由多个流域组成,每个流域具有不同的下垫面特征,因此模型参数的敏感性和最终率定值也有所不同。根据模型中各参数的敏感性大体排序,并给出参数最终率定值的范围(表2)。

模型校验效果评估表明该模型已经具备了模拟三江平原实际水文循环的能力,能够较好地表达一些难以观测的水文循环通量和状态量,为三江平原地表水地下水之间的相互作用与转化机制研究提供了有力工具。本节将利用前文构建的水文模型,分析三江平原2002—2104年期间(模型预热期不参与分析)的地表水地下水转化关系。

总体把握全区域/流域水文循环各个环节的相互转化关系有助于理解各种水文循环系统的水量平衡特点和水资源变化趋势,地表水地下水耦合模型可以较好地完成此项工作。根据水的赋存空间,将区域水文循环系统分为土壤水、地表水、地下水3大系统,各系统是相互联系、不可分割的整体,通过水分通量建立彼此的转化关系。

将三江平原地表水地下水耦合模拟输出的各项水分通量和状态变量统计到各大系统,形成补给、排泄、蓄量变化的水量平衡关系(表3),其中,土壤水系统是承接、传递水分通量最大的系统,尽管地表水系统的河道过境水量超过了2000亿m³(黑龙江、乌苏里江、松花江3大江的过境水量),但真正参与当地水文循环的仅部分人工引水、水面蒸发、渗漏等。降水是土壤水系统的主要补给项,其次是各种水源的灌溉水量,也有一部分通过毛管作用力上升至土壤层的地下水(一般称为潜水蒸发);蒸发蒸腾是土壤水系统最大的排泄项,占土壤水总排泄量的70%以上,其次是产流和深层渗漏,分别补给了地表水和地下水;从土壤水系统的蓄水量变化看,土壤层平均每年增加6.51亿m³的储量,说明降水量和灌溉量增加的趋势(图7)使三江平原的土壤趋于潮湿。

三江平原的地表水系统主要包括河流、水库、湿地等,承接的当地补给主要为产流汇入,包括地表产流、壤中流和基流;工业生活用水消耗后排放的污废水也是地表水的补给项之一;三江平原地表水的绝对损失主要是蒸发和出境,人工引水和渗漏还能继续参与当地的水文循环。总体上看,尽管三江平原建设了大量水库存蓄雨洪水,减少了出境水量,但是大规模的生产生活引水消耗导致地表水呈现逐年萎缩的趋势,各种地表水体蓄量平均每年减少近1亿m³。

三江平原地下水系统侧向相对较为封闭,与其他区域含水层基本没有侧向量交换,但垂直方向上的通量对地下水影响显著。降水向含水层的深层渗漏是地下水的主要补给方式,地表水体的渗漏补给在三江平原相对微弱,以地下水向地表水体的反向基流排泄为主。人工开采是地下水的主要排泄项,使得地下水储量以平均每年超过3亿m³的速度减少。

三大系统内部及相互间的通量交换构成区域水文循环的全过程,实际过程更加复杂,但是对于三江平原全区的 水量平衡而言,补给、排泄项相对简单。降水对于自然的封闭流域而言,是唯一的补给项。三江平原虽然是按照分水岭和三大江划分为相对封闭的区域,但又不同于封闭流域,三大江的过境水量是三江平原重要的补给项之一。三江平原通过蒸发蒸腾消耗和河道径流出境保持区域各系统水储量基本稳定。从全区域蓄水量变化上看,尽管地表水和地下水平均每年有超过4亿m³的衰减,但是总体水储量呈现增加趋势,应与三江平原降水量的增加有关。

水文模型根据下垫面、气象条件、人类活动等因素的空间差异性计算地下水数值模拟所需的源汇项通量,为平原区地下水动态演变的模拟研究提供了更加可靠的支撑。本节利用三江平原地表水地下水耦合模拟结果分析平原区地下水各项补给、排泄、蓄量变化(图7和图8)及地下水位变幅(图9)。

2002年以来,三江平原降水量呈现增加趋势。作为地下水的主要补给项,降水入渗补给量与降水量正相关。然而,作为地下水的主要排泄途径,地下水开采量也呈现显著增加的趋势。究其原因,应是降水量的增加并不能有效缓解水稻面积急剧扩张引起的灌溉需求暴涨,尽管在有些降水丰沛的年份,地下水开采量确实有所降低,但并不能改变地下水开采量迅速抬升的趋势。基流是地表水与地下水相互作用与转化的直接体现,基本上与降水量正相关,但又受到地下水开采和连枯连丰年份的影响,呈现出较为复杂的变化形势。其中,连续丰水年对基流影响较大,能使基流排泄量大幅增加。地下水开采量的持续增加导致部分地区超采严重,三江低平原中东部地区形成大片地下水降落漏斗(图9),2002年初—2014年底的13 a间,地下水降幅最大超过10 m。尽管部分地区地下水位有所上升,但并不能改变三江平原地下水储量减少的趋势。图8显示,以2000年初的地下水储量为基准,平原区地下水储量持续、大幅减少,到2008年触底,直到2011年之后储量才有所增加,但仍处于较为严重的超采状态,2014年底较2000年初储量减少46.6亿m³。

降水入渗补给系数和基流指数是表征地表水地下水相互作用与转化的2个重要系数,也是区域水资源评价中的关键参数,前者为降水量中补给地下水的比例,后者表示径流量中基流量的比例。本节将三江平原分为山丘区、三江低平原、兴凯湖平原3种特征分区,对比2个系数在不同地形地貌、水文地质、地下水开发利用程度等条件下的差异性,为不同特征区域的水资源评价、地下水数值模拟提供关键参数方面的参考。

三江平原地表水地下水耦合模型在山丘区采用松散耦合模拟方式,地下水为集总式模拟,在刻画地下水位和侧向流动方面存在不足,但在地下水相关通量模拟上,与地下水数值模拟采用的是同一套数据,都是来自子流域尺度的水文循环模拟;在对比上述2个系数时,需要把数值模拟的降水入渗补给量和基流量换算到子流域尺度上,以保证不同区域的对比具有相同的基础。

地形地貌和地质特征对降水入渗补给的影响最为显著。山丘区陡峭,降水极易形成径流汇入河道,通过岩土层渗漏补给地下水的水量相对有限,而平坦的平原区更有利于降水的滞留与下渗,提高了降水入渗补给的强度。图10a显示,三江低平原和兴凯湖平原的降水入渗补给系数远高于山丘区,前者多年平均值分别达到了0.21和0.28,而山丘区仅为0.07。即使都是平原区,三江低平原和兴凯湖平原的系数也存在差异。引起这种差异的原因可能包括2个方面：一是三江低平原中东部地层覆盖5~20 m的亚粘土,不利于土壤水的下渗,更多降水以蒸发、地表径流的方式转移,补给地下水的水量相应减少;二是兴凯湖平原区大部分年份降水量明显高于三江低平原,土壤含水量相对丰富,少量降水就能使土壤水饱和,更多的降水量下渗补给地下水。

三江平原径流量和基流指数的区域差异性极为明显,应是多种因素综合作用的结果。一般而言,平原区的径流系数(径流量与降水量的比值)远低于山丘区,如模型计算的三江低平原多年平均径流系数为0.06,山丘区为0.35,与各自的地形地貌和地质特征基本相符合。而与三江低平原自然地理特征相似的兴凯湖平原,径流系数却高达0.22,导致其径流量远大于三江低平原的径流量(图10b)。造成这种现象的可能原因,一是其年均降水量比三江低平原高13%,有更多的降水转化为了径流;二是地下水埋深相对更浅,地下水基本没有超采,降水无需弥补地下水的亏空,下渗的水量基本转化为了基流,形成径流量的主要组成部分。兴凯湖平原基流指数高达0.52,而三江低平原仅为0.24,这在一定程度上解释了兴凯湖平原径流偏多的现象。另外,山丘区基流指数始终处于较低水平,多年平均值仅为0.06,远低于平原区。由此可见,影响径流量和基流指数的因素除了地形地貌和地质特征,降水量、地下水埋深等因素也扮演着重要角色。

地表水地下水半松散耦合模型汲取了紧密耦合模型和松散耦合模型的优点,在大型流域/区域的水文循环模拟中能够较好地兼具精度与效率2个方面的优势,是目前地表水地下水耦合模拟研究的热点,也是一种较有实践应用前景的水文模型。本文介绍了一种地表水地下水半松散耦合模型的新型耦合方法,将三维有限差分地下水流数值模型嵌入半分布式水文模型中,摒弃了目前常用的HRU—网格单元空间和信息转换方式,通过建立子流域—网格单元的空间关联和地表水文—地下水流之间的实时信息交互,实现大型山丘平原交错的地表水地下水耦合模拟。

本文在三江平原构建了模型并开展应用研究。模型校验结果显示,在模拟三江平原地表水文和地下水两方面都取得了良好的效果,具备了水文循环现状模拟和未来预测的能力。本文分析了三江平原2002—2014年的地表水地下水转化关系,发现：① 土壤水、地表水、地下水之间发生着强烈的相互作用与转化,导致土壤水蓄量增加,地表水和地下水蓄量减少,但全区域总体蓄水量略有增加;② 降水和农业灌溉需求的增加推动平原区地下水的降水入渗补给和开采量整体抬升,但前者并不能填补地下水大规模超采引起的地下水储量衰减;③ 山丘区降水入渗补给系数和基流指数均显著低于平原区,但三江低平原和兴凯湖平原之间也有很大的差异,说明这2个系数不只因地形地貌和地质特征而异,还与气候气象、水文地质、地下水动态等因素有关。

模型在地表水地下水耦合模拟上取得了一定进展,今后将继续强化地表水地下水相互作用的物理性,如以数值模拟的方式刻画河流—含水层、湖泊—含水层之间的相互作用。不确定性是当前水文循环模拟面临的重要问题之一,将尝试评估影响模型的各方面不确定性,以期在模型研发和构建工作中减少不确定性方面的困扰。MODCYCLE模型与其他相似模型的对比分析是下一阶段的主要研究内容之一,通过对比明确模型在精度和效率上的优势。