入湖流量、风场等在湖体的水体运动过程中起很重要的作用，三维环境流体动力学模型EFDC提供了一个高分辨率的工具来研究太湖的水体交换.水龄可以定量反映水体的运动和交换程度以及滞留情况，对评估水质的变化具有重要意义. Deleersnijder等^[1]介绍了基于对流-扩散输运方程的水龄理论，这个新的水龄定义提供了一个健全的方法来定量地研究欧拉坐标框架下的物质溶解的输运特征. Li等^[2]运用EFDC三维水动力模型，借助水龄的定义，研究了“引江济太”引水量对太湖水循环的影响.

太湖作为中国的第三大淡水湖，由于富营养化而严重影响了饮用水的安全，政府部门近年来采取了许多措施来提升太湖水质，包括建更多的污水处理厂，禁止使用含磷清洁剂等.此外，一个重要的生态工程——“引江济太”工程在2002年启动，“引江济太”工程将更优质的长江水通过望虞河引入太湖，最终从太浦河流出，以期通过此举来加快太湖的水循环，改善太湖的水质状况^[3-5].此外，在2014年太湖北部流域内开始疏通新沟河、新孟河两条河道，太湖将通过新孟河、望虞河补水，借助走马塘、新沟河、太浦河向长江排水，形成“两进三出”的循环系统，通过让水体动起来改善水质.本文基于EFDC三维水动力模型，建立水龄模型来研究太湖水体交换速率.基于实测径流、风场数据、引水工程的调水运行数据，对太湖进行三维水动力模拟，研究风场作用下太湖水龄分布的季节性特征以及“引江济太”工程对太湖水龄分布特征的影响.

1 研究区域

太湖(30°55′40″~31°32′58″N，119°52′32″~120°36′10″E)横跨江、浙两省，西和西南侧为丘陵山地，东侧以平原和水网为主.地处亚热带，气候温和湿润，属季风气候，夏季受热带海洋气团影响，盛行东南季风，温和多雨；冬季受北方高压气团控制，盛行偏北季风，寒冷干燥.太湖是典型的风生流湖泊，且水深极浅，湖域面积几乎全年不变.太湖水域面积为2338.1 km²，平均年出湖径流量75×10⁸ m³，蓄水量44×10⁸ m³，年降水量1100~1150 mm，年平均气温16~18℃, 平均水深1.89 m，最大水深2.6 m，年平均风速3.5~5.0 m/s.为分析“引江济太”工程建设对太湖水龄的影响，环太湖共采集12个站点进行分析对比，其中1^#、5^#~7^#站点位于引水口，1^#~4^#、12^#站点位于取水口，2^#、10^#站点为水位验证站，其余站点为流量控制点(图 1).

2 数值模型 2.1 EFDC三维水动力模型 2.1.1 三维水动力模型的建立

EFDC模型是一个多参数的有限差分模型，可用于河流、湖泊、水库、湿地系统、河口和海洋等水体的水动力学和水质模拟^[6].本文将其应用于太湖的研究.控制方程在水平方向上采用笛卡尔直角坐标系或者曲线正交坐标系，垂向采用Sigma坐标.在Sigma坐标系下EFDC模型的连续性方程和动量方程为^[7-9]：

$ {\partial _t}\left( {mH} \right) + {\partial _x}({m_y}Hu) + {\partial _y}({m_x}Hv) + {\partial _z}\left( {m\omega } \right) = {Q_H} $

(1)

$ \begin{array}{*{20}{c}} {{\partial _t}\left( {mHu} \right) + {\partial _x}({m_y}Huu) + {\partial _y}({m_x}Hvu) + {\partial _z}\left( {mwu} \right) - (mf - u{\partial _y}{m_x} + v{\partial _x}{m_y})Hv = }\\ { - {m_y}H{\partial _x}\left( {g\zeta + p} \right) - {m_y}({\partial _x}h - z{\partial _x}H){\partial _z}p + {\partial _z}(m\frac{{{A_v}}}{H}{\partial _z}u) + {Q_u}} \end{array} $

(2)

$ \begin{array}{*{20}{c}} {{\partial _t}\left( {mHv} \right) + {\partial _x}({m_y}Huv) + {\partial _y}({m_x}Hvv) + {\partial _z}\left( {mwv} \right) - (mf - u{\partial _y}{m_x} + v{\partial _x}{m_y})Hu = }\\ { - {m_x}H{\partial _y}\left( {g\zeta + p} \right) - {m_x}({\partial _y}h - z{\partial _y}H){\partial _z}p + {\partial _z}(m\frac{{{A_v}}}{H}{\partial _z}v) + {Q_v}} \end{array} $

(3)

$ ({\tau _{xz}},{\tau _{yz}}) = {A_v}{H^{ - 1}}{\partial _z}\left( {u,v} \right) $

(4)

$ {\partial _z}p = - gHb = - gH(\rho - {\rho _0})\rho _0^{ - 1} $

(5)

$ ({\tau _{xz}},{\tau _{yz}}) = ({\tau _{sx}},{\tau _{sy}}) = {c_{\rm{s}}}{(U_{\rm{w}}^2 + V_{\rm{w}}^2)^{1/2}}({U_{\rm{w}}},{V_{\rm{w}}}) $

(6)

式中，x、y为水平笛卡尔坐标系；z为竖向Sigma坐标；H为水深，ζ是水面高程，η为相对水位，定义ζ=H+η；u、v、w分别为x、y、z方向的流速；Q_H为体积源项，包括降雨、蒸发、渗流；Q_u、Q_v为x、y方向的动量源项；p为附加静水压力；f为柯氏参数；τ_xz、τ_yz为x、y方向的垂向剪切应力；A_v为垂向湍流动量扩散系数；ρ和ρ₀分别为实际和参考密度，b为浮力. m_x、m_y为水平坐标因子，且m=m_xm_y. U_w和V_w为水平面10 m高度处风速x和y方向的分量，风应力系数为：

$ {c_{\rm{s}}} = 0.001{\rm{ }}\frac{{{\rho _{\rm{a}}}}}{{{\rho _{\rm{w}}}}}[0.8 + 0.065{(U_{\rm{w}}^2 + V_{\rm{w}}^2)^{1/2}}] $

(7)

式中，风速单位为m/s，ρ_a和ρ_w分别代表空气和水的密度.

2.1.2 模型主要计算参数以及模型的率定

本文采用2014年太湖的流量、气象、水位以及地形资料，基于EFDC建立三维水动力模型，EFDC模型的应用已十分成熟，时间步长通过CFL(科朗数)稳定条件确定.糙率与水深、床面形态等因素有关，根据计算率定.紊动粘性系数采用Smagorinsky公式计，C_s为Smagorinsky系数，通常取0.28^[8].数学模型的初始条件一般分热启动和冷启动两种.前者一般通过计算得到计算区域内相对准确的水位、流速场，后者则假定一种均匀分布，与实际情况不一定相符，但通过边界信息的传入，误差会逐渐消失；本文采用了冷启动.数学模型的边界条件包括开边界(水边界)和固定边界(岸边界).因入湖河流的断面尺度相对于湖域面积很小，故本文的计算将太湖视为一个封闭的区域，入湖河流的流量通过源汇项输入.模型网格采用正交曲线网格，最大网格尺度1000 m，最小网格800 m，网格数15625个，垂向分为4层，每层的厚度由湖底地形以及水面高程来决定.模型主要受径流、表层风应力的驱动，温度在模型中视为常数.模型的率定通过验证站的实测水位进行，选取犊山闸和大浦口站的水位进行模型的率定，模拟结果与实测值吻合较好，表明模型精度满足要求(图 2).由于实际情况下太湖的换水周期为300 d左右^[3]，本文采用的流量边界基本接近实际情况，因此模拟周期为365 d.

2.2 水龄的定义

水龄被定义为某个固定的水质点从入口传输到某个指定点所用的时间，通常入口处的水龄设为0.水龄可以通过系统中示踪剂的浓度分布以及一个虚拟的水龄浓度来计算.假设系统中只有一种示踪剂，那么示踪剂浓度和水龄浓度满足以下两个控制方程^[10-16]：

$ \frac{{\partial c\left( {t,\vec x} \right)}}{{\partial t}} + \nabla \left( {\vec uc\left( {t,\vec x} \right) - \vec K\nabla c\left( {t,\vec x} \right)} \right) = 0 $

(8)

$ \frac{{\partial \alpha \left( {t,\vec x} \right)}}{{\partial t}} + \nabla \left( {\vec u\alpha \left( {t,\vec x} \right) - \vec K\nabla \alpha \left( {t,\vec x} \right)} \right) = c\left( {t,\vec x} \right) $

(9)

式中，c为示踪剂浓度，$\vec u$为速度场，$\vec K$为扩散张量，t为时间，$\vec x$为空间坐标.平均水龄a可由下式计算：

$ a\left( {t,\vec x} \right) = \frac{{\int_0^\infty {\tau c{\rm{d}}\tau } }}{{\int_0^\infty {c{\rm{d}}\tau } }} = \frac{{\alpha \left( {t,x} \right)}}{{c\left( {t,x} \right)}} $

(10)

基于EFDC的源程序，增加了水龄的数值模块.为了率定风场作用下的水龄，与Li等^[2]的计算结果进行对比率定，在SE、N风向下假设单入口流量下各测站水龄计算结果拟合结果良好(表 1)，因此本文计算的风场下水龄的精度满足要求.

2.3 计算方案

本文的计算将太湖视为一个封闭的区域，太湖的环湖河道众多，但其流量均较小，因此模型在设置流量边界时，考虑上游主要入湖河流有从西苕溪水系流入的城东港、大渲河泵站、百渎口和长兜港，下游的主要出湖河流有位于东南部的太浦河，太浦河的出水量占太湖总出水量的65 %左右^[17].由于缺乏所有站点的实测资料，根据水量平衡原理，由2014年城东港、百渎口、大渲河、浯溪桥、梅梁湖泵站、望亭站的实测资料以及2014年的出湖总流量(约为104×10⁸ m³)，推算出长兜港以及太浦河的流量.太湖是典型的风生流湖泊，风场在水体循环中起很重要的作用，为研究不同风向对太湖水体交换的影响，对2014年NCEP再分析风场资料(图 3)进行分析后，得出春、夏、秋、冬的常风向分别为SE、ESE、NNE和N，采用1994-2014年的多年平均风速.在太湖湖流入口处持续释放浓度为1 mg/L的示踪剂，在不改变流量的情况下分别对四季不同风向作用下的水龄进行模拟计算，计算方案参见表 2.

3 水龄的计算以及模拟结果分析 3.1 染色剂和流场的季节性分布特征

梅梁湖口的流速实测值在0.02~0.06 m/s之间^[17]，本模型计算结果显示流速在0.03~0.07 m/s之间，与实测值吻合较好，表明模型的精度满足要求.染色剂和流场分布参见图 4，由模拟结果可知，湖流方向随风向的改变而改变，染色剂从太湖入口处放入，随着湖流方向开始扩散.春、夏偏南向季风作用下，竺山湖、西北湖区、西南湖区的染色剂浓度较大；秋、冬偏北向季风作用下，竺山湖、湖心区以及东部湖区的染色剂浓度较大.

3.2 水龄季节性分布特征

基于染色剂模型得出太湖水龄分布(图 4).结果表明，染色剂浓度较大的地方水龄较小；相反，染色剂浓度较小的地方水龄较大.湖体水龄分布统计参见表 3.春、夏、秋、冬季的湖体的平均水龄分别为220.45、235.62、217.47、211.80 d.秋、冬季的平均水龄小于春、夏季.春、夏在偏南向季风作用下，竺山湖、西北湖区以及西南湖区的水龄较小，200 d等水龄线分布范围较小；秋、冬季在偏北风作用下，竺山湖、梅梁湖、西南湖区以及湖心区的水龄较小，东部湖区、东太湖湖区的水龄较春、夏季有所减小，200 d等水龄线分布范围较大.靠近湖流入口处的地方水龄较小，远离入口的地方水龄较大；水流流向与风向一致时水龄减小，水循环加快，反之则减慢.竺山湖是湖体水龄最小的地区，贡湖则是湖体水龄最大的地区.

3.3 “引江济太”工程对水龄分布的影响

每年太湖都有低水位的“枯水期”，往年都是通过望虞河“引江济太”，确保太湖水位保持在3 m以上.为了更好地盘活太湖水体，在从长江引水进入太湖的同时，太湖流域内在2014年开始疏通新沟河、新孟河两条河道，太湖将通过新孟河、望虞河补水，借助走马塘、新沟河、太浦河向长江排水，太湖将形成“两进三出”的循环系统，通过让水体动起来，改善水质.为防止汛期太湖水位过高，新沟河将承担太湖北部区域向长江排水的功能.同时，应急情况下也能通过新沟河向太湖调水，改善梅梁湖水质.

据实测资料分析，2014年新沟河和望虞河的引水时间并不在同一区间，2014年新沟河在6-8月进行引水20 d，共引水1.1522×10⁸ m³，而望虞河2014年1-3月、11-12月共引水130 d(图 5)，引水量10.5759×10⁸ m³.据此本文对2014年望虞河引水和新沟河引水分别进行了计算对比，在引水时放入1 mg/L的染色剂，其余湖流入口采用实测流量，不放入染色剂，通过染色剂的面积可以得出进入湖体的引水水体面积，进一步分析对比从而评估引水的效益，模拟结果参见表 4和图 6.望虞河的引水量大于新沟河，因此染色剂扩散面积较大，染色剂浓度大于0.01 mg/L的面积约为新沟河的5倍；染色剂大于0.001 mg/L的面积约为太湖总面积的1/4，因此可以看出望虞河引水效果良好，表明引水量越大，引水效果越好.但是望虞河引水量约为新沟河的10倍，染色剂面积未达到10倍，是由于新沟河入口附近水体交换较快.

为进一步评估引水方案对水龄的影响，在对实际风场、流场下太湖的水龄进行计算的前提下，分别在望虞河和新沟河引水时也放入染色剂，其余湖流入口采用实测流量，均放入1 mg/L染色剂，模拟结果参见图 7~8.结果表明望虞河引水工程可以加快贡湖的水体交换，减小该地区的水龄；从新沟河引水可以加快梅梁湖地区的水体交换，改善梅梁湖的水质，但是由于引水量太小，因此对水龄的影响并不明显.

4 结论

对于湖泊的水体交换评估，水龄是一个非常重要的指标.它反映了水质点从进入湖体到离开湖体的平均时间，并受到流量、风场等因素的影响，因此可以用来评估水体的更新速率.本文建立了基于EFDC的三维水动力模型来研究太湖水龄的分布以及“引江济太”工程对其的影响.春、夏季在偏南向季风作用下，竺山湖、西北湖区以及西南湖区的水龄较小，200 d等水龄线分布范围较小；秋、冬季在偏北风作用下，竺山湖、梅梁湖、西南湖区以及湖心区的水龄较小，东部湖区、东太湖湖区的水龄较春、夏季有所减小，200 d等水龄线分布范围较大.靠近入口处的地方水龄较小，远离入口的地方水龄较大；水流流向与风向一致时，水龄减小，水循环加快，反之则减慢.望虞河引水工程能够减小贡湖湖区的水龄，加快该湖区及部分湖心区的水循环，盘活太湖的东部地区，对太湖水污染的治理具有重要意义；而2014年新沟河工程采用的引水方案可以减小梅梁湖区的水龄，加快梅梁湖区的水循环，改善梅梁湖水质，并且通过流水带动，改善入湖口直湖港、武进港的水质.在需要引水时，新孟河、新沟河、望虞河均可引水，可以组合使用，以实现最优引水方案；需要排水时太浦河、新沟河、走马塘均可排水，太浦河是主要通道.引水工程的实施对加快整个太湖水体循环做出了重要贡献.

致谢: 感谢江苏省水文水资源勘测局无锡分局的各位同仁在野外采样工作中给予的无私帮助.