生态水文双向耦合模型的研发与应用：Ⅰ模型原理与方法

曾思栋，夏 军，杜 鸿，陈向东

（1. 中国科学院重庆绿色智能技术研究院，重庆 400714；2. 长江勘测规划设计研究院，湖北 武汉 430010；3. 武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉 430072；4. 中南民族大学 资源与环境学院，湖北 武汉 430074；5. 中国水权交易所，北京 100053）

1 研究背景

水循环作为地球系统中各过程物质交换和能量转换的重要载体和纽带，与陆地生态系统过程相互作用并发生反馈［1］。一方面水循环通过降水、下渗等过程为生态系统提供水分，影响植物光合作用碳同化、CO2排放进而影响生态系统的格局［2］；另一方面陆地生态系统中植被通过气孔行为调控植物的蒸腾，植被动态变化如叶面积指数变化影响降水截留和蒸发、能量分配，对水文过程产生影响［3-4］。研究陆地水循环与生态系统的耦合是生态水文的基本内容，准确认识生态水文的耦合机制是实现流域水资源及生态管理的重要基础［5］。目前，生态水文的耦合研究已经成为水文、气象、生态等多个学科共同关注的课题。

陆地生态水文研究的主要内容体现为植被与水分、能量和物质的耦合循环机理研究［6］，构建生态水文耦合模型是认识这些耦合机理包括水、能量、碳、氮、磷等相互作用的重要基础。陆地生态水文耦合从单向耦合向双向耦合发展，单向耦合模型主要是从水文模拟角度出发，考虑植被对水文过程的单向影响，不考虑水文过程对植被动态的反馈，如DHSVM（Distributed Hydrology-Soil-Vegeta⁃tion Model）模型［7］等。随着对生态水文过程耦合认识的不断加强，双向耦合模型的理论框架被提出，如赵风华等［8］根据生态水文水、碳循环要素在植被-土壤-大气连续体（SPAC）中的运动过程，将耦合作用分为4 个尺度的关联。在双向耦合模型的构建中，由于传统水文模型对生态过程的描述不够充分，生物地球化学模型对水循环过程的描述也欠缺，因此加强分布式水文模型与生物地球化学模型的耦合可以更好地研究陆地生态系统与水循环的相互作用［9］，如Mo等［10］构建了基于能量平衡、水循环和碳氮循环的生态水文动力的VIP模型（Vegetation Interface Processes Model）；Park等［11］建立了耦合水文模型与生物地球化学模型的CHANGE 模型，并应用于西伯利亚落叶松林的水、碳通量和能量的模拟；Shen 等［12］将PAWS（Process-based Adaptive Watershed Simulator）模型与CLM 模型（Community Land Model）耦合分析不同土地利用/植被覆盖以及地下水对蒸散发及净初级生产力的影响。

具有物理机制的陆面水文模型和生物地球化学模型分别是描述陆面水文过程及生物地球化学过程的有效工具，将二者进行完全的紧密耦合对于认识陆面生态水文相互作用及其反馈机制有着重要的基础意义。目前仍然需要进一步加强基于物理机制的生态水文过程“双向”紧密耦合模型的研究，本研究将分布式时变增益水文模型（Distributed Time Variant Gain Model，DTVGM）与生物地球化学循环模型（Carnegie-Ames-Stanford Approach Carbon Nitrogen Phosphorous，CASACNP）进行双向紧密耦合构建DTVGM-CASACNP 模型，提出了耦合模型的架构、模拟方法及主要原理和方法，为陆面水文-生物地球化学过程的耦合模拟提供基础支撑。

2 生态水文双向耦合模型架构

生态水文双向耦合模型DTVGM-CASACNP 由流域水循环模块、能量平衡模块、光合作用模块、碳氮磷生物地球化学循环模块等4 大部分组成。其中流域水循环模块主要包括垂直方向大气水-植物水-地表水-土壤水-地下水等的转化以及水平方向的径流传输模拟；能量平衡主要包括冠层辐射传输、冠层及地表能量平衡以及土壤热传导中能量转化过程的模拟；光合作用模块主要包括光合有效辐射吸收、光合作用及气孔导度等碳同化相关过程的模拟，碳、氮、磷生物地球化学循环模块主要包括植物生长，碳、氮、磷元素在不同生物地球化学库包括植物库、凋落物库以及土壤库中的相互转化过程的模拟。

DTVGM-CASACNP 生态水文双向耦合模型基于DTVGM 和CASACNP 模型构建。DTVGM 模型是夏军提出并发展的通过结合空间数字化信息将集总水文模型拓展应用到流域时空变化模拟的分布式水文模型［13］，该模型将水文非线性系统理论与分布式模型框架进行有机结合，对于不同流域、水文资料信息不完整或不确定性干扰下的水文模拟具有较强的适应能力，已经在多个流域得到较好的验证［14-15］。为了描述生态水文耦合的相互作用过程，本研究对DTVGM 进行改进，将DTVGM 中基于水量平衡的两层土壤水模块改为基于Richards方程的土壤水数值模拟模型，将基于Bagrov的蒸散发模块改为土壤-植被双源蒸散发模型，增加了能量平衡模拟以及土壤热运动模块。生物地球化学模型CA⁃SACNP是由Wang等［16］和Houlton等［17］基于CASA（Carnegie-Ames-Stanford Approach）模型发展而来的，可以模拟植被动态过程以及碳、氮、磷的循环，该模型已经用于全球或区域碳、氮、磷循环的研究，并得到了较好的验证［18］。本耦合模型中植物的光合作用采用基于瞬时尺度的Farquhar［19］光合作用模型，气孔导度模型采用具有一定物理机制的Leuning［20］冠层导度模型。

生态水文双向耦合模型各模块主要运行机制如下：（1）流域水循环模块，主要计算各单元内的冠层截留、蒸散发、产流、土壤水运动、坡面汇流以及河网汇流等自然水文过程；（2）能量平衡模块，主要计算各单元冠层及土壤各项能量收支项包括短波辐射、长波辐射，潜热通量、感热通量、土壤热通量等能量通量以及冠层和地表、土壤温度变化等；（3）光合作用模块，主要计算冠层吸收的光合有效辐射、水汽及CO2气孔导度、光合作用净光合速率、初级生产力（GPP）等碳同化相关过程；（4）碳、氮、磷生物地球化学循环模块，主要计算植物自养呼吸、净初级生产力（NPP）及其分配、叶面积指数，碳、氮、磷在植物库、调落物库以及土壤库中的相互转化通量以及异养呼吸等生物地球化学过程。

DTVGM-CASACNP 以流域水循环为纽带，耦合了伴随的能量通量以及碳、氮、磷的生物地球化学过程，通过本架构实现水、能量通量、碳、氮、磷过程的紧密耦合以及反映其相互作用及反馈的实时动态关系，可以实现陆面水文-生物地球化学过程的耦合模拟，为陆地水、碳、氮、磷循环研究及水资源及生态管理提供基础支撑。图1为耦合模型的结构示意图。

图1 生态水文双向耦合模型DTVGM-CASACNP结构示意图

3 双向耦合模拟方法

陆地生态水文的双向紧密耦合模拟体现在水、能量通量、碳、氮、磷通量在每个计算时段内的相互作用及实时反馈，这些相互作用及实时反馈的描述通过关键的耦合变量实时传递实现。水、能量通量及碳、氮、磷变量实时传递主要包括：（1）能量过程为水文过程提供冠层和地表辐射传输、土壤热通量，直接影响水文过程蒸散发、土壤水运动等；水文过程为能量过程提供蒸散发、土壤含水量，直接影响能量传输过程的潜热通量、土壤热通量等；（2）能量过程为碳、氮、磷过程提供冠层和地表辐射，冠层、地表及土壤温度，直接影响植物光合作用、呼吸作用等；碳、氮、磷过程为能量过程提供植物叶面积指数等，直接影响能量的分配；（3）水文过程为碳、氮、磷过程提供土壤含水量、叶面湿润比例等，直接影响植物光合作用、植物干物质分配等；碳、氮、磷过程为水文过程提供叶面积指数、气孔导度等，直接影响蒸散发、冠层截留等。水、能量通量及碳、氮、磷通量主要变量实时传递关系与模拟流程，如图2、图3所示。

图2 耦合模型主要变量传递示意图

图3 双向耦合模拟基本方法流程图

实时动态反馈的双向耦合模拟具体算法如下：（1）将模拟时段划分为M个时段，对每个时段分别进行水、能量通量、碳、氮、磷通量的耦合模拟；（2）设定每个时段初始冠层温度、地表温度，进行冠层和地表辐射传输计算；（3）进行植物光合作用-气孔导度-蒸腾作用的耦合模拟计算；（4）计算地表产流，土壤水、热运动；（5）计算植物、凋落物和土壤中碳、氮、磷循环过程；（6）根据能量平衡迭代求冠层、地表温度，满足收敛条件结束本时段计算，不满足收敛条件回到步骤（2）直到能量平衡收敛为止；（7）进行下时段耦合模拟计算，直到所有计算时段完成。

4 模型主要原理

4.1 地表产流过程地表水产流根据时变增益非线性产流概念计算，其理论基础是降雨-径流系统的非线性关系，表征为与前期影响雨量指数（API）有关的时变增益因子［13］。计算公式如下：

式中：Rs为地表产流，mm；Ke为反映土壤平均含水量变化的参数；Pe为有效降雨，mm；g1、g2为时变增益因子。

4.2 蒸散发过程蒸散发是陆地水汽返回大气的主要形式，直接影响地表能量的收支平衡及反馈大气的潜热通量。DTVGM 原版模型中将土壤和冠层作为一个整体并通过潜在蒸散发与土壤水分胁迫来计算实际蒸散发，而陆地生态水文耦合的关键是基于植被气孔行为的蒸腾作用与光合作用的耦合。因此本研究中将蒸散发进行分层细化，分为土壤蒸发、冠层截留蒸发和植被蒸腾，并引入莫兴国等［21］对双源蒸散发模型的改进方法，计算公式如下：

4.3 土壤水运动过程土壤水运动不仅是水循环的重要环节，对于植被的蒸腾、土壤蒸发、植被生长以及碳、氮、磷的生物地球化学循环也起着关键性的作用。为细致反映土壤水的运动规律，本研究土壤水运动基于Richards方程［22］，该方程微分形式如下：

式中：θ为土壤含水率，m3/m3；ψ为土壤基质势，m；K为土壤水力传导率，mm/d；R为根系吸水项，mm；t为时间，s；z为土壤深度，m。采用van Genuchten 公式［23］描述具有高度非线性的水分特征曲线θ(ψ) 和水力传导函数K(ψ) 如下：

式中：θr、θs分别为土壤残余含水率和土壤饱和含水率，m3/m3；α、n为水分特征曲线参数。

Richards 方程的隐式差分离散格式见方程（9），模型求解采用变时间步长法，即当迭代次数大于设定次数时，将模型的时间步长变为上一次迭代的1/2，如此反复直到前后两次土壤负压满足设定精度为止。本次研究采用Celia［24］方法以减少计算过程中的水量平衡误差，用追赶法求解该方程。

式中：j为时间；m为迭代次数；z为单元中心节点的纵坐标。

本次研究上、下边界条件分别为：

4.4 植物根系吸水植被主要通过根系从土壤中吸收水分，满足其植被的生长以及代谢等生理活动和叶片蒸腾。本研究假设根系吸水R(z,t) 是蒸腾Ec(t)以及有效根分布Le(z,t)的函数，在此基础上对根系吸水加入土壤水分约束［25］。

4.5 冠层和地表辐射传输到达地表的辐射包括太阳总辐射（短波辐射直射和散射）以及天空的长波辐射，同时地表又会反射太阳短波辐射，也会发出长波辐射。对于冠层净辐射Rnc（单位为W/m2，以下同）平衡包括太阳短波辐射Rsc(↓,表示辐射传输方向向下,以下同)、冠层反射的短波辐射Rscr(↑,表示辐射传输方向向上,以下同)、天空长波辐射Rlc(↓)，冠层向天空Rlca(↑)及向地表Rlcs(↓)两个方向发射的长波辐射以及地表长波辐射Rlsc(↑)。对于地表净辐射Rns平衡包括冠层截留剩下的太阳短波辐射Rss(↓)、天空长波辐射Rls(↓)、冠层长波辐射Rlcs(↓)，地表反射的短波辐射Rssr(↑)以及地表发射的长波辐射Rlsac(↑)。

4.6 冠层和地表能量平衡冠层和地表能量平衡方程可表示为如下式：

式中：λEc、λEs分别为冠层蒸潜热通量及土壤蒸发潜热通量，W/m2；Hc、Hs分别为冠层感热通量以及土壤表面感热通量，W/m2；G为土壤热通量，W/m2；Cc和Cgs分别是冠层及地表的容积热容量，J/（m2·℃）；Tc和Tgs分别是冠层和地表温度，℃；

其中，冠层和土壤表面的感热通量［26］采用下式计算：

土壤热通量通过对两层土壤的热量估计来计算：

式中：Tgs、Ts1、Ts2分别为在地表、土壤深度d1、d2处的土壤温度，℃；k1、k2分别是土壤深度d1、d2处的土壤的热传导系数，W/（m·℃）。

4.7 土壤热运动土壤温度直接影响到能量平衡的各分量，土柱中热量传输的方程［27］如下：

式中：C为土壤体积热容量，J/（m3·℃）；T为土壤温度，℃。初始条件可以是来自站点的实测数据，或者通过气温进行估计。上边界条件为地表温度，地表温度由地表能量平衡计算中得出。下边界条件根据实际情况可以选择固定温度Tb或者零通量边界。

4.8 光合作用与气孔导度本研究采用基于Furquhar 和Collatz［28］的光合作用模型计算植被碳同化，该模型具有一定的物理机制，得到了广泛的应用。该模型中叶子的净光合速率A（umol C/（m2·s））采用RuBP 再生限制的净光合作用率Aj、羧化速率限制的净光合作用率Ac以及最大净同化速率As的最小平滑值进行计算，即将净光合速率看作是式（18）二次多项式的最小的根：

式中，Rubisco活性限制的光合速率为：

RuBP（光限制率）再生限制的光合速率为：

最大净同化速率为：

式中：Vc,max为最大羧化作用率，umol C/（m2·s）；Kc和Ko分别为CO2和O2的Michaelis-Menten 系数，umol/mol；Γ*为CO2补偿点，umol/mol；ci和oi分别为叶肉细胞间隙的CO2和O2浓度，umol/mol；J为电子传输速率，umol electrons/（m2· s）；β1、β2为 经验参数，C3 和C4 分别 指C3 植 物 和C4 植物。

蒸腾和光合作用均受到植物气孔的开关控制，本次采用Leuning 在Ball-Berry 模型［20］基础上改进的冠层导度模型，该模型基于CO2气孔导度和水汽的气孔导度成线性关系如下：

式中：gs为水汽的冠层导度，m/s；a1为经验系数；cs为叶片表面的CO2浓度，umol/mol；D为叶片表面的饱和水汽压差，kPa；Dx为表征气孔对敏感性的经验系数；g0为最小的冠层导度，m/s。

4.9 植物生长与物候模拟植物的自养呼吸包括维持呼吸Rm和生长呼吸Rg。维持呼吸［29］计算如下：

式中：cni为植物各库的碳氮比；Ci为植物各库的含碳量，umol C/m2；ri为根茎叶呼吸速率。生长呼吸通常与扣除维持呼吸的初级生产力（GPP，umol C/（m2·s））成正比，则净初级生产力（NPP，umol C/（m2·s）及其分配如下：

式中：ac,i为NPP分配到根、茎、叶的分配系数；τP,i为植物各库到凋落物库的转化速率，s-1。耦合模型中植被的物候分为4个阶段，第1阶段从发芽到稳定生长，第2阶段从稳定生长到叶子开始凋落，第3阶段从叶子凋落开始到凋落结束，第4阶段从叶子凋落结束到发芽。

4.10 碳的生物地球化学过程植物凋落后进入凋落物库，凋落物一部分被吸收到土壤有机碳库中，一部分转化为CO2。凋落物碳库的控制方程如下：

式中：bj,i为碳从植物库i到凋落物库j的分配比例；mn为土壤无机氮限制消减率；τL,j为没有被土壤无机氮限制的凋落物分解速率，s-1。

土壤有机碳分解时，一部分被其他有机碳库吸收，一部分通过异氧呼吸转化为CO2进入大气，控制方程如下：

式中：ck,j为碳从凋落物碳库到土壤碳库的分配比例；dk,kk为碳从土壤碳库kk到土壤碳库k的比例；τs,k是各碳库的分解速率。

4.11 氮的生物地球化学过程植物有机氮库分为根、茎、叶三个库，植物将从根系吸收的氮按比例分配到各库。植物衰老过程中，一部分氮被重新吸收到活的组织中，剩下的转化到凋落物库中［30］。

式中：NP,i为植物各库有机氮含量，umol N/m2；Fn,up为根系氮吸收率，umol N/（m2·s）；rn,i是重新吸收系数。

凋落物氮库的控制方程如下。

式中：nL,str为STR库的氮碳比；τL,i为各库的分解速率，s-1；mn是氮限制因子。

土壤有机氮库的控制方程如下：

式中：ck,j为氮从凋落物库j到土壤有机氮库k的比例；dk,kk为氮从土壤有机氮库k到土壤有机氮库kk的比例；τs,k为土壤有机氮各库的分解速率，s-1。

无机氮库被认为是植物根系吸氮的主要来源，其控制方程如下：

式中：Ns,min为土壤无机氮库的量，umol N/m2；Fn,dep、Fn,fert、Fn,fix、Fn,net、Fn,up、Fn,loss分别为大气氮沉降率、施肥添加率、固氮率、净矿化率、氮吸收率、氮损失率，umol N/（m2·s）。当根系吸氮量小于氮的需求量时，NPP将会减少，根系吸氮计算如下：

式中：nPmax,i、nPmin,i分别为各植物库的最大和最小氮碳比；KN,up为经验系数。

4.12 磷的生物地球化学过程磷通过大气沉降、岩石风化进入生态系统，通过淋失和侵蚀输出系统。有机磷植物各库的控制方程如下［16］：

式中：PP,i为植物各库中有机磷的量，umol P/m2；Fp,up为植被根系的磷吸收率，umol P/（m2· s）；rp,i为凋落前磷的被植物组织重新吸收系数。

凋落物氮库的控制方程如下：

式中：pL,str为STR库的磷碳比；τL,i为各库的分解速率，s-1。

土壤有机磷库的控制方程如下：

式中：ck,j为磷从凋落物库j到土壤有机磷库k的比例；dk,kk为磷从土壤有机磷库k到土壤有机磷库kk的比例；τs,k为土壤有机磷各库的分解速率，s-1；Fp,tase为磷的生物化学矿化率。

土壤无机磷库按照植被和细菌对无机磷吸收的有效程度分为3个部分，分别为不稳定磷、吸附磷和强吸附磷。以不稳定磷Ps,lab为例，其控制方程如下：

式中：Fp,dep、Fp,wea、Fp,fert、Fp,net、Fp,tase、Fp,up、Fp,loss分别代表磷沉降率、固定率、施肥添加率、净矿化率、生物化学矿化率、吸收率、损失率，umol P/（m2· s）；Spmax为吸附磷的最大量；Kplab为吸附系数。

5 结论

（1）本研究提出了生态水文双向紧密耦合的模拟方法，通过耦合的关键变量实时传递实现水、能量通量、碳、氮、磷通量在每个计算时段内的相互作用与实时反馈。

（2）本研究提出的耦合模型能够较好地刻画水、能量及碳循环等陆面水文物理、生物地球化学过程及其相互作用，为陆面水文-生物地球化学过程的耦合模拟及应用提供基础支撑。

该生态水文双向耦合模型DTVGM-CASACNP的应用案例在下篇中给出。