1. 引言

水库与其他陆地水生生态系统类似，在其生物活动的过程中伴随着有机物的分解和温室气体的释放，释放的温室气体将减损部分水电的碳汇效益 [1] [2] 。与天然湖泊相比，水位的反复波动是水库的典型特征，水位的周期性变化将导致库岸周期性淹没或露出，形成消落带。消落带面积的变化既受季节规律的影响，也受水库调度规则的驱动。水库运行过程中，库区水体中的有机质不断沉积并封存于库岸，形成高有机质含量的沉积物 [3] 。当水库消落带随水位下降暴露于大气中时，适宜的温度、湿度和有机质含量将导致消落带的温室气体排放强度增大，将成为二氧化碳等温室气体的主要排放源 [4] [5] 。由于消落带的碳源效应远高于水库蓄水前的陆地生态系统，其在碳循环中的碳源作用不容忽视 [6] [7] 。

随着对消落带碳排放现象的相关研究逐渐深入，水库碳汇属性因消落带碳排放而有所争议。苏维词 [8] 以三峡水库为例，认为水库蓄水运行后，原先适应陆生环境生长植物种群将逐步消亡，而适应淹水环境的水生植物又因消落带的季节性出露水面而难以长期生存。消落带的碳汇功能伴随植物死亡逐渐丧失，成为减损水电的碳减排效益的碳源。Marcé等 [9] 认为忽略消落带的内陆水域碳排放量计算方法存在低估，若在全球内陆水域碳排放量的基础上将消落带纳入考虑，可能会导致年均碳排放量增加0.22 Pg C，相当于总排放量的10%。Keller等 [10] 观测了全球6794个水库消落带的碳通量，认为消落带的碳排放量与碳埋藏量之比达到2.02 (1.04~4.26)，挑战了目前将水库作为净碳汇的观点。Yang等 [11] 实测了三峡水库消落带的碳排放量约为342.67~495.96 Gg/yr，抵消了消落带植被固碳量的约80%，这意味着将消落带碳排放纳入水库生态系统的碳预算对准确水库评估碳效益是必要的。消落带不同区域的碳排放强度和面积在调度过程中波动变化，在估算碳排放量时需要确定不同调度期对应的碳排放强度和消落带面积。基于碳排放因子的水库碳通量评估方法通过移用物理特性相近的水库碳排放因子等方式估算碳排放量，有利于解决实测碳通量数据缺乏的问题 [12] 。

总之，科学开展考虑消落带碳排放的生态调度可降低消落带的温室气体排放量，提升水电的碳减排效益。然而，现有针对消落带碳排放估算的相关研究主要基于全球尺度的静态数据，忽略了实际调度过程中消落带面积和碳排放强度的动态变化。本论文聚焦于水库消落带碳排放的问题，以丹江口水库与汉江中下游为研究对象，从消落带碳排放模拟和水库生态调度模型构建两个方面开展研究，提出了基于消落带面积和碳排放强度季节差异的碳排放量计算方法，建立了考虑消落带碳排放的多目标优化调度模型，通过多目标优化求解得到了兼顾生态效益与社会经济效益优化调度图。

2. 水库生态调度模型构建

2.1. 目标函数

2.1.1. 消落带碳减排目标

1) 消落带面积计算

水库消落带在调度过程中受到水位波动影响呈现周期性的淹没或露出，其面积决定了当前时段消落带的碳排放量。利用水库的水位库容曲线，结合消落带坡度计算不同水位对应的消落带的面积 [11] 。假设相邻水位之间的库容横截面为直角梯形，相邻水位对应的库容差除以水位差得到水面面积，将最大水面面积与各水位对应水面面积作差，除以消落带坡度的余弦值即可得到该水位对应的消落带面积 [13] 。

2) 消落带碳排放强度计算

消落带碳排放强度随着淹没露出呈现周期性的变化，不同阶段有着显著的差异。根据淹水时间和范围将消落带划分为落干区、退水区和淹水区 [3] [14] 。

消落带落干区指长期未被水淹没或淹没后经历足够长时间排干的区域，伴随季节性的植被生长和演替，其碳排放强度可根据土壤呼吸模型估算 [2] [4] [15] 。土壤呼吸模型的计算参数包括土壤温度(ST)、土壤湿度(SM)、植被叶面积指数(LAI)和土壤有机碳含量(SOC)，计算公式为 [16] ：

$Rs=R_0\times {\text{e}}^{\left(\frac{\ln Q_{10}\times ST}{10}+0.093\times SM\right)}\times \frac{LAI}{LAI+0.496}\times \frac{1.876\times SOC}{SOC+1.726}$(1)

式中：Rs表示土壤呼吸强度，单位µmol/(m²·s)；R₀表示0℃时的土壤基础呼吸量，由碳–质量–温度理论确定，单位µmol/(m²·s) [17] [18] ；Q₁₀表示土壤呼吸强度对温度的敏感性，即土壤温度每上升10℃土壤呼吸强度的增加量；ST表示土壤温度，单位℃；SM表示土壤湿度；LAI为植被叶面积指数；SOC表示土壤有机碳含量，单位kg C/m²。

消落带退水区指水位下降后重新露出水面的区域。消落带淹没后，植被固存水下，水体中的有机物不断沉积，土壤有机质含量升高 [1] [19] 。同时，长期淹没后土壤的含水量高，适宜微生物的生存繁殖，导致其呼吸强度升高 [20] 。因此，退水区相较落干区碳排放强度更高，根据实测数据，水库消落区的碳排放强度约为落干区的3.78倍 [3] 。随着退水过程推移，土壤含水量逐渐降低，土壤有机质逐渐分解或转化为其他形式，土壤呼吸强度降低至落干区水平。

消落带淹没后形成淹水区，气体交换由土–气界面转变为水–气界面，土壤呼吸受到抑制。淹水区的碳排放强度主要受到水生动植物的呼吸、光合作用影响，不同季节具有显著差异 [21] 。

水库消落带的碳排放量由落干区、退水区和淹水区的面积与对应的碳排放强度决定。在调度的过程中，尽可能低的消落带年均碳排放量(AACE)意味着更好的水库碳减排效益。构建描述消落带碳减排的目标函数为：

$\text{min}AACE={\displaystyle \underset{t=1}{\overset{T}{\sum }}\left(F_t^{\text{D}}{A}_t^{\text{D}}+F_t^{\text{R}}{A}_t^{\text{R}}+F_t^{\text{I}}{A}_t^{\text{I}}\right)\Delta t/N}$(2)

式中：AACE表示年平均碳排放量，T表示总调度时段数，N表示调度年数， $F_t^{\text{D}}$、 $F_t^{\text{R}}$、 $F_t^{\text{I}}$分别表示第t时段落干区、退水区和淹水区的碳通量， $A_t^{\text{D}}$、 $A_t^{\text{R}}$、 $A_t^{\text{I}}$分别表示第t时段落干区、露出区和淹水区的面积， $\Delta t$表示各调度时段的时间间隔。

2.1.2. 鱼类适宜生境目标

1) 水文变异度最小

天然流量过程包含刺激生物节律的特征信号，为鱼类繁殖创造了适宜的环境。水库建成后，河流由激流状态转变为受人工调控的稳定流态，天然河道被渠化江段取代。通过计算不同调度方案的水文变异度，降低水库调控对天然流量过程的改变有利于鱼类种群的发展，构建水文变异度(Deviation of variability range, DVR)的目标函数为：

$\min DVR=\sqrt{{\displaystyle \underset{i=1}{\overset{I}{\sum }}\left(\frac{Y_{o,i}-Y_e}{Y_{o,i}}\right)/I}}\times 100\%$(3)

式中：DVR表示总体水文变异度，Y_o,i表示水文参数i的观测值落入RVA范围的年数；Y_e表示落入RVA范围的年数期望值(取总年数的一半)；I表示水文参数的个数。

2) 适宜鱼卵漂流的河段长度最大

除对天然水文情势的需求外，鱼卵孵化过程中的流速需求易被忽视。粘沉性卵对低流速的需求集中于2月到4月，一般要求流速不高于0.8 m/s；漂流性卵的需求时段集中于5月到8月，一般要求流速不低于0.3 m/s；其他时段鱼类不处于产卵期，可不考虑对流速的需求。考虑到仙桃产卵场受损及汉江下游鱼类资源萎缩的问题，建立一维水动力模型确定流量、流速关系，计算满足鱼类繁殖需求的河段长度。调度过程中，满足鱼类需求的流程越长，鱼卵在漂流过程中顺利孵化的可能性就越大，构建描述适宜鱼卵漂流的河段长度(SRL)目标函数为：

$\max SRL={\displaystyle \underset{t=1}{\overset{T}{\sum }}R{L}_t/T}$(4)

式中：SRL表示年均适宜鱼卵漂流的河段长度，RL_t表示第t时段的适宜鱼卵漂流的河长，T表示调度时段数。

2.1.3. 河流水华防控目标

1) 浮游植物密度最小

对汉江中下游历史水华监测资料的分析揭示了浮游植物密度和水体营养状态是水华现象的表征因素，反映了水华的生消过程 [22] 。目前针对水华防控效果的评价主要基于监测手段，缺乏完善的机理模型用于构建生态目标。因此，可利用监测资料和随机森林(RF)建立数据驱动模型，探究水华与流量等水文要素间的内在联系 [23] [24] [25] 。在水华频现期1月~3月，尽可能低的平均浮游植物密度(APD)意味着水华暴发的危害降低 [26] 。构建浮游植物密度的目标函数为：

$\min APD={\displaystyle \underset{t=1}{\overset{T}{\sum }}P{D}_t/T}$(5)

式中：APD表示研究时段的平均浮游植物密度，PD_t表示第t时段的浮游植物密度。

2) 营养状态偏离率最小

水体富营养化是水华暴发的促发条件之一，可采用综合营养指数评价(I_TL)。合理的泄水过程能够减轻河流富营养化程度，降低富营养时段的比例 [26] 。构建营养状态偏离率(DRI)目标函数为：

$\min DRI=\#\left(I_{TL,t}\ge 50\right)/T$(6)

式中：DRI表示营养状态偏离率， $I_{TL,t}$表示第t时段的水体综合营养指数。

2.1.4. 社会经济目标

1) 年均供水量最大

构建年均供水量(AAWS)最大的目标函数为 [27] [28] ：

$\text{max}AAWS={\displaystyle \underset{t=1}{\overset{T}{\sum }}W{S}_t/N}$(7)

式中：AAWS表示年平均供水量，WS_t表示第t时段的供水量，N表示调度期年数。

2) 年均发电量最大

构建年均发电量(AAHG)最大的目标函数为 [29] [30] [31] ：

$\text{max}AAHG={\displaystyle \underset{t=1}{\overset{T}{\sum }}\left(K{Q}_{\text{O},t}H\right)\Delta t/N}$(8)

式中：AAHG表示年平均水力发电量，K表示水电站出力系数，Q_O,t表示第t时段水库下泄流量，H表示发电水头。

2.2. 边界条件

水库调度过程中的水量平衡约束、水库库容约束、出库流量约束、最大出力限制约束、最小生态基流约束如下 [32] [33] ：

$V_{t+1}=V_t+\left(Q_{\text{I},t}-Q_{\text{O},t}\right)\Delta t$(9)

$\underset{\_}{V}\le V_t\le \bar{V}$(10)

$\underset{\_}{Q_{\text{O}}}\le Q_{\text{O},t}\le \overline{Q_{\text{O}}}$(11)

$P_t\le \bar{P}$(12)

$Q_{\text{e}}\le Q_{\text{O},t}$(13)

式中：V_t表示第t时段的库容，Q_I,t表示第t时段的入库流量，Q_O,t表示第t时段的出库流量， $\underset{\_}{V},\bar{V}$分别表示库容上下界， $\underset{\_}{Q_{\text{O}}},\overline{Q_{\text{O}}}$分别表示出库流量上下界，P_t为第t时段的水电站出力， $\bar{P}$为水电站最大出力限制，Q_e表示下游河道内最小生态流量。

2.3. 优化算法

布谷鸟搜索算法适用于各种优化问题，已经成功地用于解决许多工程实例。对丹江口水库3条供水调度线的拐点进行优化。为了兼顾水库调度图的准确性和计算效率，参考常规调度图，各供水调度线设置4个拐点。优化变量的上边界在汛期为防洪限制水位，在非汛期为正常高水位；优化变量的下边界为死水位。优化的3条供水调度线不应相互交叉。设置多目标优化模型的种群规模为100，迭代次数为2000。

3. 研究对象

3.1. 汉江中下游流域

汉江中下游干流总长约652 km，自西北向东南纵贯湖北省28个县市区，是湖北省经济发展的核心地区之一 [34] 。汉江中下游的七级梯级水利枢纽建成后，兴隆枢纽以上河段成为梯级水库相连的渠化河段，部分指示水生生物生命节律的重要特征信号消失，出现了鱼类生态群落萎缩、渔获物小型化趋势明显、鱼类产卵场破坏的问题 [35] [36] [37] 。汉江中下游流量呈明显季节趋势，在冬末春初进入枯水期，硅藻等浮游植物易在气象水文条件适宜的条件下暴发式增长，发生水华，威胁河流生态健康 [22] [38] 。

3.2. 丹江口水库

丹江口水库(110˚59'E~111˚49'E; 32˚36'N~33˚48'N)位于汉江中游，水面横跨湖北省丹江口市和河南省南阳市，由1973年丹江口大坝下闸蓄水后形成，具有供水、防洪、水力发电、灌溉、航运等综合效益。丹江口水库水质优良，被作为南水北调中线工程的水源地，承担了向河南、河北、北京、天津等四省市的供水任务，设计年均供水量达95亿m³[39] 。丹江口水库供水调度图的时间尺度为旬，有五个运行区 [40] [41] 。丹江口水库属于平原型水库，库岸坡度较缓。坝顶加高后水库淹没面积达到1153.7 km²，在调度的过程中会产生大面积的消落带，其温室气体排放会减损部分水力发电的碳减排效益 [1] 。

4. 生态调度结果

4.1. 消落带碳排放目标构建

4.1.1. 丹江口水库消落带面积变化

丹江口水库属于平原型水库，库岸坡度平缓，具有较大的消落带面积。根据叶松等 [42] 的研究结果，丹江口水库消落带150~160 m平均坡度约为4.21˚，160~170 m平均坡度约为5.97˚。采用水位库容曲线估算水库消落带的面积，结果如图1所示。

图1. 水位库容曲线和消落带面积水位关系

丹江口水库消落带面积与水位呈负相关，在死水位145 m时最大，达552.00 km²。采用水库库容曲线关系估算的消落带面积与基于DEM模型的计算结果基本一致，差距仅为0.14% [42] [43] 。

4.1.2. 碳通量强度及季节分布

淹水区的碳通量强度由设置在丹江口水库水面的通量站实测，落干区的碳排放强度基于碳排放因子估算得到，年内变化如图2所示。可以看出，丹江口水库水面碳排放强度具有明显的季节趋势，夏季碳源，冬季碳汇，年内碳源和碳汇效应基本抵消，年均碳排放强度为6.93 mmol/(m²·d)。消落带落干区和退水区全年作为碳源，碳排放强度随温度呈现明显的季节性变化，年均碳排放强度约为195.68 mmol/(m²·d)。落干区的土壤呼吸强度夏季较高，在8月达到峰值，约323.06 mmol/(m²·d)。退水区的土壤呼吸强度根据落干区估算，在退水一段时间后逐渐降低，并于30日后降低至落干区水平。

为了验证采用本论文方法估算消落带碳排放强度的合理性，与以往研究成果进行比较。Yang等 [11] 曾实地监测水库消落期的土–气界面CO₂通量，约为188.09 mmol/(m²·d)；Li等 [7] 通过逐月采样和野外测量，测得消落带的土–气界面的CO₂通量约为249.12 mmol/(m²·d)，与本论文成果基本一致。

图2. 消落带各区域碳排放强度年内分布

4.2. 多目标优化调度

4.2.1. 优化参数设置

本论文采用布谷鸟多目标算法对丹江口水库3条供水调度线的拐点进行优化。为了兼顾水库调度图的准确性和计算效率，参考常规调度图，各供水调度线设置4个拐点。优化变量的上边界在汛期为防洪限制水位，在非汛期为正常高水位；优化变量的下边界为死水位。优化的3条供水调度线不应相互交叉。设置多目标优化模型的种群规模为100，迭代次数为2000。

4.2.2. 优化调度方案

根据目标函数和边界条件建立优化调度模型，制定考虑消落带碳排放的丹江口水库生态调度方案。通过多目标优化调度，得到一系列帕累托解。图3展示了帕累托解集的平行坐标图。

图3. 帕累托解集的平行坐标图

对比帕累托解集与常规调度结果，筛选各目标均不劣于常规调度且总提升率最大的方案为优化方案。图4展示了与优化调度方案相对应的水库调度图，按照优化调度图进行水库调度，可得出对应的供水流量、下泄流量、水库水位的变化过程和各目标值。

表1展示了优化调度与常规调度的结果，以及各目标的改进度。优化调度方案能够兼顾生态效益与社会经济效益，在保证AAWS和AAHG不低于常规调度的前提下，改善生态目标，总改进度达8.41%。AAWS和AAHG改进度达1.13%和1.86%，依照优化调度图开展调度能够提高水库的社会经济效益，提高水资源利用率。AACE

图4. 优化调度图

表1. 常规调度方案与优化调度方案结果对比

的改进度最大，达3.19%。优化调度方案调整了不同时段的供水任务，从而改变了消落带碳排放量的季节分布，可见调度方案差异对消落带碳排放量的影响较为显著。DVA降低了1.04%，SRL提高了0.52%，这意味着优化调度方案的泄流过程更符合流域天然水文情势，适宜鱼类种群的繁衍和恢复。尽管DVA和SRL的提升程度不是太大，但其小幅度提升对鱼类而言仍至关重要，这是由于适宜的短期流量事件就能够促使鱼类捕捉到指示其生命周期的水文信号。ADP和DRI都有所降低，但是改进度较小。出现该现象的原因是针对浮游植物密度的研究时段集中于河流枯水期，此时能够用于调度的水量较少，不同调度方案的出库流量过程都接近最小生态流量，优化空间较小；优化调度图未改变水库泄水的季节规律，各月的综合营养指数主要受流量大小的影响，改进度较小。

4.2.3. 调度目标间的协同竞争关系

各目标之间的相关性反映了社会经济目标和生态目标之间的竞争协同关系，可为水库调度方案的实施提供参考。计算优化方案之间的相关系数r，结果如图5所示。

图5. 各调度目标间的相关性

如图5所示，AACE与AAWS、ADP呈协同关系，与AAHG、SRL、DRI呈竞争关系，与DVA相关性较弱。为探究各目标在不同供水调度方案下的结果，以年均供水量为基准绘制其他目标结果的散点分布，如图6所示。

图6. 以AAWS为基准的调度目标散点图

丹江口水库采取库区引水的方式向南水北调中线工程供水，供水量无法被水电站再次利用，AAWS与AAHG呈竞争关系 [32] 。当AAWS增加时，AAHG的损失也逐渐加快，这是由于供水量增大导致水库水位下降，上下游水位差降低，在减少相同水量的情况下发电效益受损更严重 [27] [30] 。

AACE与AAWS呈协同关系，这主要包括两方面的原因：丹江口水库库区引水直接减少水库蓄水量，扩大消落带面积；出库流量(包括供水流量和下泄流量)与入库流量的不匹配会造成水位波动，使得退水区面积增加，消落带碳排放量增加。值得注意的是，随着AAWS增加，AACE的增速加快。当AAWS小于90.72亿m³时，AACE增速缓慢，此时AAWS与AACE之间的竞争较小，落干区和退水区面积主要随来水的季节变化波动；当AAWS进一步增大，供水任务与AACE之间的竞争加强，AACE随AAWS的增长速度加快。因此，在满足受水区需求的前提下，选取供水量较低的调度方案将提高水库调度的碳效益 [12] 。

DVA呈现随AAWS的增加先减小后增大的趋势。当供水量较小时，IHA子指标落入RVA范围的比例高于期望值，导致DVA较高；随着供水量逐渐增大，各子指标落入RVA范围的比例低于期望值，DVA呈上升趋势。因此，保持中等强度的供水量有利于降低DVA，恢复适宜鱼类种群的流量特征。SRL与AAWS呈现负相关关系。水库供水量增加意味着生态水量的减少，若河段在孵化期未达到适宜流速阈值将导致漂流性鱼卵下沉死亡，影响鱼类种群的维持 [35] 。

APD随AAWS增加呈现上升趋势。由于浮游植物密度随季节、流量的变化较为剧烈，不同调度方案的泄水时机不同，对浮游植物密度影响较大，导致数据点较为散乱。而随着AAWS的增加，在水华暴发期可用于调度水量减少，水库维持最小生态流量下泄，流量过程平缓，无法抑制浮游植物的繁殖。DRI随AAWS的增加而降低，即加大水库泄流将提高营养状态破坏率。高流量导致的水体富营养化主要集中于夏季，此时江汉平原的农业活动频繁，河流淹没范围加大，裹挟土壤有机质、化肥等面源污染物输入河流，水体理化因子较高。

5. 结论

本论文针对水库消落带的碳排放问题，以丹江口水库与汉江中下游为研究对象，从消落带碳排放模拟和水库生态调度模型构建两个方面开展研究。论文主要研究内容与结论如下：首先，针对水库消落带的碳排放问题，提出了考虑消落带面积和碳排放强度季节差异的碳排放量计算方法；其次，建立了考虑消落带碳排放的多目标优化调度模型，求解得到帕累托解集，提出了兼顾经济效益和生态效益的水库优化调度图；最后，根据调度目标之间的相关性明确了其协同竞争关系。研究结果表明：基于消落带面积和碳排放强度的碳排放量计算方法反映了消落带碳排放的季节趋势。丹江口水库消落带面积与水位之间的相关性强，消落带落干区和退水区全年作为碳源，淹水区夏季作为碳源，冬季作为碳汇，年内的碳源效应和碳汇效应基本抵消。优化调度方案能够兼顾生态效益与社会经济效益，在保证AAWS和AAHG不低于常规调度的前提下，改善各生态目标，其中年均消落带碳排放量降低3.19%。AACE与AAWS、ADP呈协同关系；与AAHG、SRL、DRI呈竞争关系；与DVA相关性较弱。本论文提出的模型和方法适用于不同的水库调度场景，有助于指导水库调度规则的制定，为水资源管理者提供参考依据。

基金项目

嘉陵江亭子口水库以下梯级联合生态调度项目(LZ-DL-2022-050)。

参考文献