气候变化和人类活动对鄱阳湖水龄影响的定量区分*

范宏翔，徐力刚,3，朱 华,2，鲁 照,4，曹宇贤,2，吴亚坤,4，姜加虎

(1：中国科学院南京地理与湖泊研究所，中国科学院流域地理学重点实验室，南京 210008)(2：中国科学院大学，北京 100049)(3：中国长江三峡集团有限公司长江生态环境工程研究中心，北京 100038)(4：安徽工业大学能源与环境学院，马鞍山 243000)

数值模拟被认为是预测水体中污染物的迁移转化的重要手段之一[2-3].截至目前，已有大量研究依托数值模型在鄱阳湖开展水环境的研究工作(表1)，如李冰[4]基于MIKE21构建鄱阳湖二维物理过程模型，研究了水情条件变化对湖泊防洪、水资源、水环境，乃至湖泊总体水安全的影响，结果表明影响鄱阳湖水安全状况的主要因素是流域来水变化、长江干流来水变化、流域入湖污染物通量和湖区人类活动.最近，国内部分学者针对设想中的鄱阳湖水利枢纽工程对鄱阳湖湖区水环境的可能影响进行了模拟预测研究，结果表明在水利枢纽不同运用方案下，枯水期主湖区的总氮和总磷浓度低于现状，但尾闾区域污染物浓度显著升高，加之湖流较现状条件下变缓，鄱阳湖局部地区的富营养化风险增加[5].综上，大多数研究主要着眼于鄱阳湖水环境本身的变化及外部其影响因素，而对于其内在机制的研究较少.特别是在气候变化和人类活动共同作用下，湖泊内部的水动力条件也发生了巨大的变化，其对水环境的影响机制尚不清楚.

表1 鄱阳湖数值模拟研究进展

为了定量认识流场在物质输运中的作用，有必要研究与流场密切相关又能反映其本质的辅助变量，比如水龄，进而通过研究这些辅助变量，有助于人们更好地理解湖泊水动力过程的环境效应.水龄作为换水能力的一种重要表达，可以定量反映湖泊水体的交换能力和滞留状况，对评估水质的变化以及富营养化风险识别具有重要意义[10, 18].鄱阳湖的水龄研究较少，并且也大多都集中在水龄对水环境的影响方面，如黄爱平[10]就基于水龄提出了鄱阳湖富营养区域的识别方法；范宏翔[19]则发现水文情势对鄱阳湖水质的影响主要通过影响湖泊的水龄来体现.但是，对于水龄本身变化及其影响因素的研究较少.气候变化和人类活动影响着全球和区域水文循环过程，是导致水文水资源时空分布不均的重要因素，同时也是流域-湖泊水文情势变化的主要原因[20].定量区分人类活动与气候变化对湖区水动力条件的影响分量，对于湖区水环境污染风险识和预警别具有重要的科学意义.目前，已有大量研究区分了二者对流域水文情势的影响[21-23]，然而对于气候变化和人类活动对湖区内部水动力条件的影响研究则十分缺乏，这主要是由于其对水文数据有着较高的要求，需要不同情景下长时间尺度连续的水位流量数据.近期，深度学习网络的发展解决了时间序列的长短期依赖问题，大大提高了时间序列预测的精度.其在水文领域的应用也越来越广泛[24-27]，利用深度学习网络构建流域径流预测模型，并通过与历史时期的数据进行对比，可以得到精细的时间尺度上的径流变化数据，并为进一步驱动水动力模型提供数据基础.

综上所述，本文目的旨在(1)构建鄱阳二维水动力模型与气象径流模型，模拟计算鄱阳湖水龄的时空变化过程；(2)通过引入基准期概念，并耦合构建的二维水动力模型与气象径流模型，利用差异化的情景模拟，定量区分鄱阳湖水龄改变的主要影响因素.研究结果能够揭示鄱阳湖水环境对水文水动力变异的响应，为鄱阳湖水污染防治提供科学支持，同时为流域水资源管理提供一种新的思路和方法.

1 材料与方法

1.1 基础数据

根据资料搜集情况，选取鄱阳湖流域时间序列最长的14个国家级气象站(图1)的日气象数据(包括降雨、温度、气压、湿度、日照时数、风速、饱和水汽压等)，数据时段为1960-2015年.鄱阳湖入湖水量主要来自“五河”，分别为赣江、抚河、信江、饶河和修水.选取5个子流域下游出口断面流量监测站(图1)日平均径流数据来训练和验证鄱阳湖流域降雨-径流模型.

图1 鄱阳湖流域地理位置及气象站分布

1.2 研究方法

1.2.1 鄱阳湖流域气象-径流模型构建 基于长短记忆模型框架构建神经网络模型来模拟鄱阳湖流域的径流过程.构建的模型有一个输入层，共包含128个LSTM神经元，用来输入流域的气象参数.在此模型中，输入层为经过标准化后的鄱阳湖流域14个国家气象站点的逐日观测数据，包括日降雨、日平均气温、日最高气温、日最低气温、日平均气压、太阳辐射、相对湿度、极大风速以及近地面风速.另外，为了防止模型出现过拟合的现象，在输入层之后设置了一个丢弃层，并设置其丢弃率为0.4[28].模型输出层是一个包含5个独立输出端的全连接层.输出数据为经过标准化后的赣江、抚河、信江、饶河和修水的逐日径流过程.模型主要的参数类型包括全局超参数、训练时超参数以及普通参数等，其中，全局超参数以及训练时超参数通过经验值进行设定，普通参数在模型训练时进行自动优化.模型结构及主要参数设定参见范宏翔等[29].

1.2.2 鄱阳湖湖区水动力模型构建 本文采用EFDC(environmental fluid dynamics code)模型构建鄱阳湖湖区二维水动力模型.EFDC模型是美国环境保护局推荐使用的水环境生态模型之一[30]，最早由美国弗吉尼亚海洋研究所开发[31]，其物理机制和很多计算方法借鉴了曾被广泛使用的Blumberg-Mellor模型[32]和美国陆军工程兵团开发的切萨皮克湾模型[33].EFDC模型包括了水动力、风浪、泥沙、富营养化、有毒物质和沉积物多个模块，并良好地耦合为一个整体[30]，可以进行水体一维、二维和三维流场、盐分、黏性和非黏性泥沙输运、生态过程和淡水入流的模拟.模型在垂直方向采用σ坐标变换，水平方向采用直角或者曲线正交坐标系.水动力模块是EFDC模型的核心，其水动力方程采用有限差分法进行求解，水平方向采用交错网格离散，时间采用二阶精度的有限差分法和内外模式分裂法进行积分.其中，外模块采用半隐式格式，允许较大的时间步长；内模块使用垂直扩散的隐式格式.EFDC在潮间带区域采用干湿网格技术，能够很好的刻画水面的干湿变化过程.EFDC模型应用广泛，能够模拟河流、河口、湖泊、水库、湿地和自近岸到陆架的海域，并能同时考虑风、浪、潮、径流和水工建筑物的影响.

1.2.3 模型基本设定和概化 在模型构建的基本设定和概化思想中，应充分考虑以下因素：(1)鄱阳湖湖底地形起伏较大，湖中主河道形状弯曲且深度较大，水动力条件复杂，网格构建时应该对主河道区域有足够的刻划；(2)鄱阳湖属于亚热带湿润季风气候区，有较为明显的季风特征，春季和冬季受到西伯利亚冷气流的影响，盛行偏北风；夏季和秋季受到太平洋副热带高压控制，盛行偏南风，多年平均风速为3.01 m/s.尽管风作为影响鄱阳湖水文情势变化的次要因素，对整体水位模拟效果和精度影响较小，但是风场也是湖泊水动力过程的重要驱动力，在某些时刻对湖泊水流结构会产生重要影响[34].因此，在构建鄱阳湖二维水动力模型时，必然要考虑风场条件对鄱阳湖水动力模拟过程的影响；(3)由于鄱阳湖“高水是湖，低水似河”的独特水文特性，在构建鄱阳湖二维水动力耦合模型时，认为水动力指标在垂向近似均匀，以垂向平均值表达[8].鄱阳湖流域气象径流模型和湖区二维水动力模型耦合过程如图2所示.

根据美国采暖、制冷与空调工程师学会（ASHRAE）提出的标准，“冷水机房全年综合能效”在0.85kW/ton以下的为高效能机房，即COP值4.4以上，在1.0kW/ton以下的为需要改造的机房，即COP值3.5以下。可见本机房能效高于国内平均水平，在低负荷工况下，其COP值将更高。冷冻机系统运行能耗情况如表3所示。

图2 鄱阳湖气象径流模型和水动力模型耦合框架

1.2.4 模型地形网格剖分 模型采用2010年鄱阳湖湖盆地形原始资料(基面为85国家高程)，通过导入基于测量和数字化的地形数据来设置湖盆地形.模型计算域边界则根据鄱阳湖淹没范围确定，在水平方向将计算区域或划分为16736个计算单元，采用曲线-正交网格，单元尺度为：i方向(东西向)约56～1578 m，j方向(南北向)约30～1420 m.由于鄱阳湖湖区地形起伏较大，为了更准确地模拟鄱阳湖的水动力过程，对鄱阳湖主河道区域的网格进行了局部加密.主河道区域的网格分辨率大致为30 m×60 m(图3).

图3 模型所用地形和计算网格

1.2.5 模型初始条件与边界条件 模型选取2012年全年水位、水质过程进行率定验证.为了减少初始条件对模型的影响，尤其是初始流速场对模拟结果的影响，在本研究中，模型采用热启动的方式进行启动，即启动前先预热1个月，水动力和水质模块的预热期为2011年12月1日至12月31日.正式模拟时，水动力模型的初始水位场和流速场采用预热后的水位场及流速场作为输入.鄱阳湖的主要入湖河流共有5条(简称“五河”).本模型共定义了8个开边界，其中包括7个流量边界和1个水位边界.7个流量边界作为水动力模型的上游开边界，分别为赣江主支、赣江中支、赣江南支、抚河、信江、饶河和修水，其表征了流域来水对湖泊水文水动力水环境的影响.将北部湖泊出口-湖口的水位过程线作为模型的下游开边界，用来表征鄱阳湖与长江之间的出流、顶托和倒灌效应.上游各个开边界的流量过程根据五河尾闾的主要控制站(外洲、李家渡、梅港、渡峰坑、虎山、虬津和万家埠)和各河道从控制站到入湖的传播时间给定[4].其中，赣江各条支流流量根据其占外洲站流量的份额按经验分配.此外，由于鄱阳湖有25000 km2的区域属于无站点控制的集水区，为保证湖泊水量平衡，提高水动力模块的准确性和合理性，需要合理估算未控区间的产流过程，并与流域“五河”控制站的流量过程共同作为模型的上游边界条件[4, 15].本文采用降雨-径流系数法对其进行估算：

W=P·ΔA·a

(1)

式中，W表示某年未控区间的产流总量(m3)；P表示该年降雨总量(m3)；ΔA表示鄱阳湖流域未控区间的集水区面积(km2)；a表示鄱阳湖流域的评价径流系数, 在本文中取经验值0.6[35].

1.2.6 模型评价指标 模型的模拟效果采用常规的统计指标来量化评价，在本研究中主要采用的指标为纳什效率系数(NSE)、均方根误差(RMSE)以及平均绝对误差(MAE)，其计算公式分别为：

(2)

(3)

(4)

1.2.7 模型数据可视化 基于R语言中的ncdf4、data.table以及tidyverse软件包，独立开发了EFDC模型前后处理软件包efdcr(https://hxfan1227.github.io/efdcr/).efdcr在极大程度解决了EFDC模型结果的可视化问题，具有极强的灵活性.同时，efdcr还提供了一系列工具函数来辅助进行EFDC的模型构建过程.efdcr的主要功能和工作流程如图4所示.

图4 efdcr工作流程

1.2.8 情景设置方案 为了定量区分气候变化和人类活动对鄱阳湖水龄变化的贡献率，本文设置了3个情景.S0作为基准情景，即自然状况下鄱阳湖流域的径流过程，在设置时参考前人的研究，将1960-1969年作为基准期[23,36-37].同时，考虑到LSTM模型的特殊性，基准期选取时也考虑了不同的流域来水情景.根据水文情报预报规范(GB/T 22482-2008)，按照鄱阳湖流域年最大流量进行水文年的划分.其中，1962年为典型丰水年，1963年为典型枯水年，1966年和1968年为典型平水年.因此，我们认为选用1960-1969作为基准期，能够满足情景模拟的需要.为了进一步区分人类活动与气候变化对鄱阳湖径流的影响，在模型训练和验证时仅使用基准期的气象水文数据，为了防止数据泄露问题，将基准期划分为训练期和验证期两个阶段：模型的训练周期从1960年1月1日至1968年12月31日，共计3287 d.模型的验证周期从1969年1月1日至1969年12月31日，共计365 d.S1鄱阳湖气象-径流模型模拟的1970-2015年鄱阳湖流域“五河”多年平均日流量过程，S2为实测的1970-2015年鄱阳湖流域“五河”多年平均日流量过程.通过对比S1和S0，可以定量区分气候变化对鄱阳湖水龄变化的影响分量；而通过对比S2和S1，则可以定量区分人类活动对鄱阳湖水龄变化的影响贡献.

1.2.9 影响程度计算 参照Ye等[23]关于定量区分气候变化和人类活动对鄱阳湖流域径流量影响程度的分析方法，由气候变化和人类活动引起的鄱阳湖水龄变化可以用下式来解释：

(5)

式中，WA0、WA1和WA2分别代表S0、S1以及S2情景下鄱阳湖水龄的变化过程，ΔWAC和ΔWAH分别代表气候变化和人类活动引起的鄱阳湖水龄变化的绝对分量.

2 研究结果

2.1 鄱阳湖气象径流模型验证

训练期模型在各个子流域的NSE分别为赣江0.94、抚河0.95、信江0.95、饶河0.95和修水0.94；RMSE分别为赣江630.80 cm3/s、抚河143.64 cm3/s、信江194.74 cm3/s、饶河129.46 cm3/s和修水144.01 cm3/s.而在验证期，模型在各个子流域的NSE分别为赣江0.90、抚河0.95、信江0.95、饶河0.98和修水0.96；RMSE分别为赣江671.03 cm3/s、抚河158.90 cm3/s、信江212.04 cm3/s、饶河 139.52 cm3/s和修水162.79 cm3/s(表2).因此可以认为构建的气象-径流模型能够很好地模拟鄱阳湖流域的径流过程.

表2 鄱阳湖流域气象-径流模型模拟表现

2.2 鄱阳湖二维水动力模型验证

验证选取了鄱阳湖湖区4个水文站点(自北向南依次为星子、都昌、棠荫和康山).验证时段选择平水年2012年1月1日至12月31日.其中，水位验证的结果如图5所示.如图所示，本文构建的鄱阳湖水环境模型的水动力模块模拟效果良好，其模拟的水位序列与实测序列吻合度较高.

图5 鄱阳湖湖区各水位站水位验证结果

由鄱阳湖2012年水位模拟精度可以发现，4个水文站点水位模拟的决定系数(R2)均在0.97以上，均方根误差(RMSE)在0.30～0.58 m之间，纳什效率系数(NSE)均在0.96以上.平均绝对误差(MAE)在0.20～0.50 m之间(表3).综上所述，本文构建的鄱阳湖水动力模块能够较好地反映鄱阳湖湖区水量交换及平衡的季节性变化特征，模拟结果真实可靠.

表3 2012年鄱阳湖水位模拟精度

2.3 基准情景下鄱阳湖水龄的时空分布特征

自然条件下(S0)鄱阳湖水龄的时空分布特征如图 6所示.可以发现，鄱阳湖湖区的水龄存在明显的空间分异特征，其最主要特征为东部湖区和南部湖区尾闾水龄较大，平均为228.01 d，而大部分通江水体水龄较小，平均为24.21 d.在洪泛湿地区域，水龄有比较明显的时间差异.枯水期洪泛湿地与主湖断开连接，形成独立的水体，导致其水龄增大，年平均水龄约为90 d左右.而在丰水期，随着水位上涨，洪泛湿地逐渐与主湖连接，形成一个整体，其水龄逐渐减小，平均为30 d左右.

图6 自然条件下(S0)鄱阳湖水龄时空分布特征

2.4 气候变化与人类活动对鄱阳湖水龄的影响分析

图 7和表4为气候变化对鄱阳湖水龄影响程度的时空差异(S1～S0)和影响范围统计，由图7可以发现，在气候变化的影响下，鄱阳湖全湖大部分湖区的水龄都呈现下降的趋势，下降幅度较小，约为20%左右；而在湖区东部，有部分区域水龄呈现上升的趋势，最大上升幅度超过100%.

图7 气候变化对鄱阳湖水龄影响程度的时空差异

从影响范围来看(表4)，在气候变化影响下，年均有2054 km2的水体水龄呈现减少的趋势，但是减少幅度较小；相反，年均有736 km2的水体水龄呈现上升的趋势，上升幅度约在30 d左右.

表4 气候变化对鄱阳湖水龄影响范围统计(km2)

进一步区分了人类活动对鄱阳湖水环境的影响程度.人类活动的影响计算方式为(S2-S1)/(S2-S0).图8和表5为人类活动对鄱阳湖水龄影响的时空差异和影响范围统计.由图8可以发现，人类活动对鄱阳湖水龄的影响存在明显的时空差异.从空间上来说，在人类活动影响下，东部湖区的水龄呈现显著的下降趋势，最大下降幅度超过30%，而中、北部湖区水体水龄则总体呈现上升的趋势，最大上升幅度超过100%.其余部分的水体水龄变化并没有呈现明显的时空特征，变化并无规律.

图8 人类活动对鄱阳湖水龄影响程度的时空差异

表5 人类活动对鄱阳湖水龄影响范围统计(km2)

从影响范围来看(表5)，在涨水期和丰水期(3-7月)仅有545 km2的水体水龄呈现下降趋势，占全湖水面积的13%，下降程度在0～32 d之间，同时，有超过2372 km2的水体水龄呈现上升趋势，约占全湖面积的58%.而在枯水期和退水期，平均有718 km2的水体水龄呈现上升趋势，占全湖面积的18%.

3 讨论

本文耦合了深度学习网络和传统二维水动力模型，通过引入基准期概念，定量区分了气候变化和人类活动对鄱阳湖湖区水龄变化的贡献程度.在以往的研究中，流域径流的模拟往往依托传统的水文模型，如SWAT[38]和WATLAC[13]等.但是，流域模型构建的过程中，所需要的数据资料庞杂，且往往涉及环保、水文、气象、农林各个部门[38]，难以全部获取，限制了物理模型在数据缺乏地区的使用[28, 40].而本文提出利用流域气象数据结合深度学习网络模拟流域径流过程，同时与湖泊二维水动力模型进行耦合，也能很好地反映流域水文状况对湖泊的影响，同时，其大大降低了模型对流域资料的依赖性.近年来，随着很多高级的神经网络框架如Keras和Pytorch被开发出来，构建一个复杂的神经网络正在变得越来越方便，为在更加精细的尺度上区分气候变化和人类活动对水文变量影响提供了数据支持和方法借鉴.

尽管深度学习网络在模拟时间序列方面具有很大的优势，模拟精度也相对较高，但是仍有大量学者认为此类黑箱模型因为不具备物理机制，无法对结果进行解释[41-43].最近，Kratzert等[44]通过对比LSTM中神经元的状态，可以揭示某些气象变量的控制机制，但是合理性还有待于验证.如何科学解释神经网络的结果也成为了学界研究的重点.在本研究中，我们发现鄱阳湖水龄主要受气候变化的影响，其相对贡献率超过了80%.气候变化通过改变流域产流过程，进而影响湖泊水龄变化.有研究指出流域入流是影响鄱阳湖水龄的主要因素，流域入流量越大则鄱阳湖水龄越小.当流域入流增加20%时，湖口水龄均值减少1.15 d[10].在本研究中，气候变化使得鄱阳湖水龄均呈现下降的趋势，其根本原因可能是气候变化使得鄱阳湖流域的降雨增加，从而导致流域入湖径流量提高.同时，蒸发作为影响流域水量平衡的另一重要因素，其变化对流域径流的改变也有着决定性作用.有研究显示，近年来鄱阳湖流域的潜在蒸发呈现长期下降趋势[23]，而实际蒸发也呈现出波动下降的趋势[45].蒸发的下降也会在一定程度上导致鄱阳湖流域径流的增加.综上，在气候变化作用下，鄱阳湖流域入湖径流量大量增加，导致其水龄呈现下降的趋势.

本研究同时也区分了人类活动对鄱阳湖水龄变化的影响分量，可以发现，人类活动对鄱阳湖水龄变化也有显著的影响.人类活动通过影响流域来水，从而间接影响了鄱阳湖湖区的水龄变化.一般来说，人类活动对流域水文过程的影响主要表现在两个方面，即改变流域下垫面特征以及兴建大型水利工程[46].以江西省为例，有研究表明随着江西山江湖开发治理工程的实施[47]，从1985-2001年，全省水土流失面积从330万hm2下降到130万hm2；全省森林覆盖率由31.5%上升到59.7%.植被覆盖的大量增加，使得流域径流迅速减少，径流的减少会导致鄱阳湖湖区水龄呈现一定的上升趋势.同时，我们发现人类活动对水龄变化的影响存在一定的年内时空差异，主要体现为枯水期水龄增大，丰水期水龄减小.这可能是由于水利工程的修建导致流域径流变化呈现年内分布不均的趋势.有研究指出[48]，三峡水库调度对中下游水情影响主要表现为11月至次年4月径流量加大，5-10月径流量减小.3月径流量增加最大，10月减少最多，9月也减少较为明显，这很好地揭示了人类活动对水龄影响的年内差异.但是不可否认的是，本文没有考虑流域下垫面因素对径流的影响.尽管本文参考前人研究选取1960-1969年作为基准期，来尽可能减少人类活动的影响，但是有研究显示，鄱阳湖流域自1950s开始就经受了剧烈人类活动，从而也对流域水循环造成了一定的影响[23].同时，鄱阳湖流域修筑有大量的水库，这些水库的蓄水也会对流域的产流过程造成一定的影响[49].在未来的研究中，应该综合考虑这些要素，将流域下垫面情况进行参数化建模，使得模型能够更加准确地反映气候变化和人类活动对流域水文过程的影响.

4 结论

1)本文基于EFDC构建的鄱阳湖二维水动力模型，能够较好地反映鄱阳湖湖区水量交换及平衡的季节性变化特征，模拟结果真实可靠.4个水文站点水位模拟的R2均在0.97以上，NSE均在0.96以上.

2)鄱阳湖湖区的水龄存在明显的空间分异特征，其最主要特征为东部湖区和南部湖区尾闾水龄较大，平均为228.01 d，而大部分通江水体水龄较小，平均为24.21 d.在时间尺度上，鄱阳湖水龄的差异主要体现在洪泛湿地区域.枯水期洪泛湿地年平均水龄约为90 d左右.而在丰水期，洪泛湿地年平均水龄为30 d左右.

3)通过耦合鄱阳湖流域气象-径流模型与湖区水动力模型，定量区分了人类活动与气候变化对鄱阳湖水龄变化的贡献率.在气候变化的影响下，鄱阳湖年均有2054 km2的水体水龄呈现减少的趋势，减小幅度为30 d左右；人类活动对水龄变化的影响存在一定的年内时空差异，在涨水期和丰水期有约58%的湖区水龄呈现上升趋势，而在枯水期和退水期，平均有718 km2的水体水龄呈现上升趋势，占全湖面积的18%.从影响权重来看，气候变化是造成鄱阳湖水龄变化的主要因素，其影响权重为84%，而人类活动的影响仅占16%.