鄱阳湖湖泊流域系统水文水动力联合模拟*

李云良，张 奇，姚 静，李相虎

(1:中国科学院南京地理与湖泊研究所湖泊与环境国家重点实验室，南京 210008)

(2:中国科学院大学，北京 100049)

鄱阳湖位于长江中下游南岸，在优质淡水供给、洪水调控和生物多样性保护等方面起着重要作用[1].鄱阳湖流域主要由赣江、抚河、信江、饶河和修水五大子流域和湖区平原区构成.鄱阳湖接纳流域来水，调节后由北端湖口注入长江(图1a).流域内降雨和蒸发具有显著的时空变化，流域人类活动强烈，极端气候事件频发，多种因素的相互叠加，致使流域内洪涝灾害频发[2-4].特别是极端气候事件的发生，影响了流域水资源的时空分布和湖泊水量变化.已有研究表明，流域五大水系的入湖径流量对湖泊水位起着主控作用，其次，鄱阳湖与长江季节性的水量交换也对湖泊水位有着不同程度的影响作用[5-6].湖泊水位的变化，导致鄱阳湖水面面积呈现明显的萎缩与扩张，严重威胁湖区平原区农业生产和管理及湿地生态系统健康等方面[5].鄱阳湖水情和水位变化，主要受季节性降雨和流域人类活动及长江来水等多因素影响，这些因素的叠加影响了流域降雨-径流过程和湖泊泄水条件，从而影响湖泊的水量平衡[6-7].由此可知，鄱阳湖湖泊流域之间具有显著的水文水动力联系，主要表征为流域来水对鄱阳湖水文的驱动作用.在目前极端气候事件频发和人类活动加速背景下，研究鄱阳湖水位和水量对流域径流过程的响应成为一个紧迫的任务.为此，本文提出一个联合模拟方法，将作为一有效工具模拟分析鄱阳湖水文水动力对流域径流的响应.

鄱阳湖湖泊流域系统自然属性复杂，如流域尺度较大(16.2×104km2)，多河入湖，湖泊岸线及湖盆地形复杂，湖泊水动力过程受流域与长江的共同作用，这些都为水文水动力过程的模拟带来挑战.水文或水动力模型往往在单一的流域或地表水体(水库、湖泊等)具有较强的模拟能力，但多数情况下与周围水系割裂，无法切实模拟复杂水文水动力系统[8].本文所提出的基于不同模拟功能的子模型的有效联合拟解决鄱阳湖湖泊流域的水文水动力整体模拟.针对鄱阳湖及其流域，已有众多研究成果且强调不同的目的[6-7，9-11].大部分研究主要采用水文模型侧重于流域降雨-径流过程的模拟计算[7，10]及其对气候变化和人类活动的响应[6]，但均未涉及到湖泊的相关研究以及湖泊对流域来水的响应，由于流域与湖泊作用关系密切，单方面的流域模拟研究显然是不足的.叶许春采用分布式水文模型WATLAC 探讨未来不同的气候变化情景对鄱阳湖流域径流的可能影响，继而采用湖泊水量平衡原理，计算流域径流变化对湖泊水位的影响.该方法没有考虑湖泊的水动力特性，只能在平均意义上描述湖泊水位的变化，存在一定的局限性[11].国内外关于鄱阳湖湖泊流域系统的水文水动力联合模拟研究几乎没有.本文的主要研究目的是，通过不同子模型的联合来实现大尺度鄱阳湖湖泊流域系统的水文水动力整体模拟;阐述该联合模拟模型的驱动方法、率定及其潜在的用途.本文提出的方法也可用于其它类似湖泊流域的水文水动力过程研究，为定量模拟湖泊水情对气候变化和人类活动的响应提供有效工具.

1 材料与方法

1.1 基础数据与准备

模型需要的主要数据包括:流域水文站点河道日径流量数据、湖泊站点日水位数据、日气象资料(降雨和蒸发皿数据)、土地利用数据、土壤物理属性数据(田间持水量、饱和渗透系数及总孔隙度)、叶面指数(来自MCD15A2)、地面高程数据(DEM)、分辨率为30 m×30 m 湖泊地形高程数据以及流域与湖泊边界矢量图等资料.降雨和蒸发站点(图1a)数据按最近邻法插值到流域和平原区计算单元，作为水文和平原区产流模拟的驱动因素;水文模型中土地利用(2005年)划分6 种类型:耕地、林地、草地、水体、居民用地和其它非建筑用地;土壤总孔隙度、田间持水量、饱和渗透系数等均为1 km×1 km 空间栅格数据，可直接输入模型，作为土壤物理属性参与计算;基于遥感反演的叶面指数数据，因其在日尺度上变化相对较小，水文模型按月尺度输入计算.湖泊水动力模型所需的主要数据为湖泊地形数据和岸线边界(图1b)，水动力模型同样按最近邻法获取站点的降雨和蒸发资料.

1.2 模型描述

图1 鄱阳湖湖泊流域系统结构(a:湖泊流域平原区结构简图及水文气象站点空间分布;b:鄱阳湖岸线边界及湖泊地形)Fig.1 Lake Poyang-catchment system (a:schematic diagram for lake-catchment-plain area and spatial distribution of hydro-climate stations;b:shorelines and bathymetry for Lake Poyang)

WATLAC 模型[12-13]是基于格网的以日为步长、联合地表-地下径流的分布式流域模拟模型，能针对湖泊集水域的特点，将相对独立的多个子流域、多条入湖河流在同一个模型中加以模拟，无需对每个子流域分别建立独立的模型，并便于对整个流域进行全局模型参数率定.其模拟的主要水文循环过程为:冠层截留与蒸散发、坡面流、河道汇流演算、土壤壤中流、土壤蒸发与渗漏补给、地下水蒸发、河流与地下水的相互作用等.有限差分方法MODFLOW-2005[14]的植入使得该模型能够刻画较为复杂的地下水流环境，是一个物理机制较为明确的地表水-地下水流实时耦合的流域模拟模型.详细计算方法这里不再阐述，具体的模型耦合机制及数学方程已于 2009年发表[12-13]且成功应用于不同的流域[7，9-13].

平原区产流的计算主要基于传统的水文学方法，采用基于概念性的平原区产流估算方法，即降雨乘以径流系数(流域多年统计结果)来获取产流量.平原区这种近似估算方法虽难以得到有效验证，但该方法已成功应用于其它流域[15]，在缺乏资料的情况下，可作为平原区入湖径流的近似估算.本文主要目的并不是关注于该方法计算的平原区入流过程的精度，至少提出的概念性评估方法弥补了以往平原区产流计算的空白，确保联合模型的完整性及结构的合理性.

MIKE 21 模型[16]是丹麦水力研究所(DHI)开发的系列水力学软件之一，属于平面二维水流数学模型，该模型可用于模拟河流、湖泊、河口、海湾、海岸等大型水体流场特征.因其具有完全的物理机制且具有灵活的三角形网格剖分，所以能够很好适应鄱阳湖的宽浅型地形特点和复杂岸线特征.该模型已成功应于不同国家及地区[17-18]，详细的模型理论及方法，这里不再阐述.

为了保证模型空间离散分辨率与大多数基础数据初始分辨率一致且能更好反映下垫面自然属性特征，WATLAC 模型的离散空间分辨率选择较为精细的1 km×1 km 格网单元，共剖分138634 个计算单元，计算域(图1a，DEM 覆盖区域)占整个湖泊流域系统面积比重约为87%.模型中仅考虑潜水含水层的模拟且概化为一层，主要用来计算河流与含水层的水量交换.其它模型详细的构建过程与Ye 等方法相同[9]，主要不同之处除了当前流域模型采用更为精细的格网外，最主要的是当前模型充分考虑平原区这部分产流的模拟计算.

平原区产流模拟计算域(图1a，灰色区域)同样采用1 km×1 km 格网，共剖分17556 个单元，占整个湖泊流域系统面积比重约11%.由于鄱阳湖流域降雨-径流关系比较密切，本文采用降雨乘以多年平均径流系数(α=0.58)的方法来估算平原区入湖径流量.径流量的估算大致分为以下两步:首先，将日气象站点的降雨资料按最近邻法插值到整个平原区格网单元，将格网单元的降雨量乘以多年平均径流系数，便可得格网的产流量;其次，假定平原区沿着湖区边界的分布式入流与点入流有着等同的影响效果，故将计算格网产流采用最近邻原则分配至5 个子流域出口，逐时段计算，便可得五大子流域与平原区的日合成径流序列过程，充分体现全流域模拟的完整性.

鄱阳湖水动力模型计算域范围约为东西向100 km，南北向170 km，占湖泊流域系统的面积比重约2%.模型采用基于有限体积法的三角形网格以适应复杂的湖泊地形及岸线边界(图1b).灵活地变网格大小为70 ～1500 m，共剖分网格数为20450.网格的合理剖分及大小不仅适应湖区内河道深窄且流速较快、平原宽浅且水流较慢的地形特点，也提高了水动力模型的计算效率.为了更加合理的模拟，模型将站点降雨、蒸发资料考虑到水动力模型中参与计算.以湖泊各站点水位空间插值结果作为初始水位场，初始流速设定为0.为了合理描述湖泊地形特点，变曼宁数给定为30 ～50 m1/3/s 分别适应于湖泊底部的平原区和主河道.涡粘系数采用Smagorinsky 公式[19]计算，其中Smagorinsky 因子Cs 取为0.28.模型中通过定义临界水深来确定干湿单元，很好地去适应湖泊水面面积随水位变化显著的特点.

1.3 模型框架

为了充分体现大尺度复杂联合系统的完整模拟，湖泊水动力模型的上游边界选在5 大子流域(赣江、抚河、信江、饶河、修水)出口断面处(图1b)，以5 大子流域的日模拟径流(约占湖泊88%贡献量)及平原区入湖径流(约占湖泊12% 贡献量)的合成流量以5 个点入流的方式作为水动力模型的上边界条件，其反映了流域-湖泊相互作用关系;模型下边界采用实测湖口水位过程线，其反映了长江-湖泊相互作用关系;其它均按闭边界条件给定.3 个子模型的连接采用输入-输出的外部耦合技术将流域WATLAC 模型、平原区产流模型以及湖泊水动力模型MIKE 21 联合起来，形成一个流域-平原区-湖泊结构的水文水动力联合模拟框架(图2).

图2 湖泊流域系统水文水动力联合模型驱动关系及模拟框架Fig.2 Driving relationship and framework of integrated hydrological and hydrodynamic models for Lake Poyang-catchment system

2 率定策略

选用2000年1月1日-2005年12月31日作为整个联合模型的率定期.通过6 个水文站点(图1a)的河道日径流量来率定WATLAC 模型的地表水流部分.将上述6 个水文站点的日径流量进行基流分割[20]，并采用加权平均值法计算整个流域的均基流指数(基流量/河道流量)，以此来评估地下水模拟效果.总之，以6个站点的河道径流量和地下水基流指数作为多目标函数来率定分布式水文模型.

在当前联合模型参数率定方面，PEST[21]参数自动率定技术主要用于WATLAC 模型，这是因为其能够与WATLAC 模型较为方便地建立数据及参数传递接口.PEST 是一个非线性参数估计工具，其原理是采用GML算法来优化模型，通过尽量少的迭代和模型调用次数使得计算值与观测值之间的残差平方和达到最小，以此来寻求最优参数值.PEST 自动率定技术的嵌入使用，能够提高模型的运行效率和整体性能，加速模型率定进程和增强模拟结果可靠性，同时也为WATLAC 模型大尺度流域、长时间连续模拟提供保障.WATLAC模型率定的主要敏感性参数为坡面流滞后系数、马斯京根法的时间蓄量常数、地下水补给率、降雨入渗因子及潜水含水层渗透系数和给水度[12-13].鉴于PEST 优化技术对模型参数的率定已取得成功应用[10]，故本文没有列出具体的参数优化结果.而水动力模型MIKE 21 仍然采用手工试错法率定，模型可调整的主要参数为湖床糙率系数和涡粘系数.

3 结果与讨论

3.1 模型率定与评估

当前联合模型采用基于观测数据来独立率定各个子模型.水文水动力联合模型分别以五河6 个水文站点河道径流量、流域基流指数与湖泊4 个站点水位作为目标变量来评价模型的整体能力.

表1 为基于多目标变量的WATLAC 模型定量化模拟结果评价.6 个水文站点拟合的纳希效率系数[22](Ens)变化范围为0.71 ～0.84，确定性系数(R2)介于0.70 ～0.88 之间，相对误差(RE)基本控制在±10%内，但个别站点拟合的相对误差稍微偏大(表1).模拟值与观测值大部分集中在45°线附近，表明模型模拟效果较好(图3).此外，通过图4 选取的模拟期末2005年河道日径流量过程线可见，总体呈现较好的曲线拟合宏观效果，但个别时段洪峰流量拟合存在一定的误差，这主要是因为流域存在众多大小型水库改变了径流的天然状态.WATLAC 模拟的全流域基流指数(45%)与基流分割结果(47.6%)较为接近，表明地下水对河道径流基流贡献量的模拟结果具有一定的可靠性，基流指数可用来直接率定分布式水文模型及间接评估地下水流模拟效果.总体来说，流域水文模型WATLAC 能够很好地再现大尺度流域上的降雨-径流响应过程，能为湖泊水动力模型提供可靠的流量输入条件.

表1 流域水文模型WATLAC 率定结果*Tab.1 Calibration results of catchment hydrological model WATLAC

图3 模型率定期(2000-2005年)模拟与观测的河道日径流量散点图Fig.3 Scatter diagrams of simulated stream flow against the observed stream flow for calibration period 2000-2005

湖泊MIKE 21 模型的模拟与观测水位的定量化评估结果表明，湖泊各站点拟合的确定性系数R2介于0.96 ～0.98，纳希效率系数Ens变化范围为0.88 ～0.98，相对误差RE 均小于±5%，取得令人满意的率定结果(表2).鄱阳湖4 个水文站点实测与计算水位过程线拟合良好，均很好地呈现了湖泊水位的季节性变化和年际变化趋势(图5).总体来说，湖泊水动力模型MIKE 21 不但能够取得理想的模拟效果，而且具有较强的模拟能力去再现长时间的水位变化及捕捉水位的峰值变化，充分体现该模型能够很好的模拟湖泊水位对流域和平原区径流量的共同响应过程.尽管如此，在枯水位的模拟上有着相对较大的误差(图5)，误差来源最可能是湖盆地形的误差或者说地形误差对低水位的模拟较为敏感.

图4 2005年鄱阳湖流域6 个水文站点模拟与观测的河道日径流量拟合Fig.4 Comparison of simulated and observed daily stream flows at six hydrological stations in Lake Poyang catchment in 2005

图5 2000-2005年鄱阳湖4 个水文站点水位验证Fig.5 Validation of lake levels at four hydrological stations in Lake Poyang during 2000-2005

图6 2005 年4 个典型时段模拟鄱阳湖水位的空间变化Fig.6 Simulated spatial distribution of Lake Poyang water levels at four time slices in 2005

表2 湖泊水动力模型MIKE 21 率定结果*Tab.2 Calibration results of lake hydrodynamic model MIKE 21

图7 2005年鄱阳湖淹水天数统计结果Fig.7 Statistical result of submerged days of the entire Lake Poyang area in 2005

3.2 鄱阳湖空间水位模拟

鄱阳湖水位在年内的空间上存在着明显的水位梯度，即水位空间差异性较大，特别是低水位月份，比如1月份(图6a)，只有湖泊主河道有少量的水分布，大部分区域水深较浅，呈现出露状态.而水位相对较高的7月份(图6c)，整个湖泊保持着较高的水位，且整个湖区近似呈水平面.此外，本文将2005年的逐日空间水位分布结果进行统计，将水深大于零的格网单元定义为湖区淹水状态，反之，定义为出露状态，空间结果见图7.鄱阳湖整个湖区的淹水情况同样呈现较大的差异性，不同的湖区有着不同的淹水情况(图7).常年被水淹没的区域主要分布在主河道地区，该淹水分布与低水位的空间分布相似(图6a)，表明这些主河道区域常年有水的流动.而湖泊岸线附近的上游区域则呈现1 ～6 个月不等的淹水时间，表明这些地区由于湖盆地形较高，一年中的绝大部分时间则处于出露状态.图6 和图7 充分体现了鄱阳湖在不同季节河湖相交替的独特水情动态，也表明当前所构建的联合模拟模型能够合理呈现鄱阳湖水位在年内的空间变化，而这种水位空间动态变化及淹水信息对评估湖泊湿地系统、水生植被的生长和发展等有着重要的现实意义.

4 结语

鄱阳湖湖泊流域系统涉及尺度大、多个独立子流域入湖、湖泊水动力特征显著等特点，系统内水文水动力过程复杂.本文尝试了一个实现该系统水文水动力的联合模拟的方法.该方法基于不同功能子模型的有机联系，通过数据在模型间的时空传递实现了湖泊和流域的联合模拟.该联合模拟方法在鄱阳湖湖泊流域作了验证，能有效反映湖泊水位对流域径流过程的响应.本文的研究虽然针对鄱阳湖湖泊流域系统，所提出的方法具有一定的普适性，可为相似湖泊流域的水文水动力模拟提供借鉴或被直接采用.

[1]周文斌，万金堡，姜加虎.鄱阳湖江湖水位变化对其生态系统影响.北京:科学出版社，2011.

[2]Shankman D，Qiao LL.Landscape changes and increasing flood frequency in China's Lake Poyang region.The Professional Geographer，2003，55(4):434-445.

[3]Shankman D，Heim BD，Song J.Flood frequency in China's Lake Poyang region:trends and teleconnections.International Journal of Climatology，2006，26:1255-1266.

[4]Shankman D.River management，landuse change，and future flood risk in China's Lake Poyang region.River Basin Management，2009，7(4):423-431.

[5]Hu Q，Feng S，Guo H et al.Interactions of the Yangtze River flow and hydrologic processes of the Lake Poyang，China.Journal of Hydrology，2007，347:90-100.

[6]郭 华.气候变化及土地覆被变化对鄱阳湖流域径流的影响[学位论文].南京:中国科学院南京地理与湖泊研究所，2007.

[7]刘 健，张 奇，左海军等.鄱阳湖流域径流模型.湖泊科学，2009，21(4):570-578.

[8]Xu ZY，Adil NG，Thomas JG.The hydrological calibration and validation of a complexly-linked watershed-reservoir model for the Occoquan watershed，Virginia.Journal of Hydrology，2007，345(3/4):167-183.

[9]Ye XC，Zhang Q，Bai L et al.A modeling study of catchment discharge to Lake Poyang under future climate in China.Quaternary International，2011，244:221-229.

[10]Li XH，Zhang Q，Xu CY.Suitability of the TRMM satellite rainfalls in driving a distributed hydrological model for water balance computations in Xinjiang catchment，Lake Poyang basin.Journal of Hydrology，2012，426/427:28-38.

[11]叶许春.近50年鄱阳湖水量变化机制与未来变化趋势预估[学位论文].南京:中国科学院南京地理与湖泊研究所，2010.

[12]Zhang Q，Li LJ.Development and application of an integrated surface runoff and groundwater flow model for a catchment of Lake Taihu catchment，China.Quaternary International，2009，208(1/2):102-108.

[13]Zhang Q，Werner AD.Integrated surface-subsurface modeling of Fuxianhu Lake catchment，Southwest China.Water Resour Manage，2009，23:2189-2204.

[14]Harbaugh AW.MODFLOW-2005:The U.S.Geological survey modular ground-water model-the ground-water flow process.Geological Survey Techniques and Methods，2005.

[15]Kebede S，Travi Y，Alemayehu T et al.Water balance of Lake Tana and its sensitivity to fluctuations in rainfall，Blue Nile basin，Ethiopia.Journal of Hydrology，2006，316:233-247.

[16]DHI.MIKE 21 ＆ MIKE 3 FLOW MODEL FM Hydrodynamic Module.Water ＆ Environment ＆ Health，2007.

[17]Niemann SL，Jensen JH，Zyserman JA et al.Morphological modeling of a danish tidal inlet.Proceedings of ICCE，2006:1-14.

[18]Martinelli L，Zanuttigh B，Corbau C.Assessment of coastal flooding hazard along the Emilia Romagna littoral，IT.Coastal Engineering，2010，57:1042-1058.

[19]Smagorinsky J.General circulation experiment with the primitive equations.Monthly Weather Review，1963，91(3):99-164.

[20]Arnold JG，Allen PM，Muttiah R et al.Automated base flow separation and recession analysis techniques.Ground Water，1995，33(6):1010-1018.

[21]Doherty J.PEST:Model-independent parameter estimation.Watermark Numerical Computing，2005.

[22]Nash JE，Sutcliffe IV.River flow forcasting through conceptual models part 1—a discusstion of principles.Journal of Hydrology，1970，10:282-290.