牛栏江滇池补水工程是滇池流域水环境综合治理六大工程措施的关键性工程，是滇中调水的近期重点工程，可有效增加滇池水资源总量和提高水环境容量，加快湖泊水体循环和交换。自2013年底工程正式通水以来，截至目前已累计向滇池补水约23亿m³，有效改善了滇池水质^[1]。该工程主要包含德泽水库水源枢纽工程、干河提水泵站工程及输水线路3部分。利用干河泵站从德泽水库取水后，水体经输水线路到达位于盘龙江松华坝水库下游2.2 km处的瀑布公园，然后经盘龙江汇入滇池。但由于德泽水库的入库泥沙偏细，而取水口以上的水库回水距离较短，加之取水口上游河段频繁发生的异重流现象及干、支流交汇作用对含沙量垂向分布的扰动等，使得实际运行中存在汛期取水含沙量明显偏高的问题，导致进入瀑布公园、盘龙江和滇池的水体明显浑浊，对滇池水环境及瀑布公园和盘龙江水景观带来了较大影响^[1-2]。

针对牛栏江滇池补水工程汛期存在的泥沙问题，前期在深入了解影响牛栏江滇池补水工程汛期取水含沙量因素的基础上，结合水库泥沙及取水防沙治理的经验^[3-4]，本文提出了降低取水口含沙量的新建分层取水口和修建拦沙坝2种方案^[1-2]。新建分层取水口是通过改造当前取水口来完成，而修建拦沙坝则是在取水口上游的干、支流内分别修建拦沙坝。

取水口布置是决定取水工程成功与否的关键因素，当取水水质要求较高，且水体的温度、浊度、盐度等随水深出现较大变化时，应考虑分层取水方式，根据取水要求抽取不同高程的水体^[5]。因此，为了满足供水水质、灌溉水温及水库下游水生态发展的要求^[6-8]，分层取水得到了广泛的应用。一般而言，分层取水可利用含沙量的垂线分布规律，取到表层较清的水体^[9]，但对水沙运动较为复杂的取水河段，分层取水方式对降低取水含沙量的效果并不明确。拦沙坝的作用主要是通过壅高水位，降低水流流速，将入库泥沙拦截在库内，从而达到减少泥沙进入下游河道的目的。以往关于拦沙坝的研究多针对沟道中拦截泥石流及对水土保持的作用^[10-13]，而对于水库中拦沙坝拦沙效果对下游取水口作用的研究^[14]、特别是对拦截细颗粒泥沙效果的研究较少。

SCHISM三维水沙数学模型^[15-16]能够满足河流—河口—海洋跨尺度计算的要求，适应剧烈变化的地形条件，在维持较高计算精度的同时采用高性能并行计算的程序架构，大幅度提高计算效率。该模型已成功运用于美国切萨皮克湾^[17]、波罗的海和北海连接区^[18]、葡萄牙海岸系统^[19]等沿海地区和黄河下游河道^[20]的计算，具有较好的适用性和准确性。考虑到德泽水库地形复杂及入库泥沙较细的特点，同时为便于更好地描述德泽水库的异重流现象，本文采用SCHISM模型研究不同典型年的汛期水沙过程在德泽水库中的运动规律，对比分析修建拦沙坝前、后沿程及取水口分层含沙量的变化情况，探讨修建拦沙坝以及新建分层取水口方案对降低取水口含沙量的影响效果。研究成果既可为制定牛栏江滇池补水工程汛期泥沙问题的合理治理方案提供参考，又对研究细颗粒泥沙河道中修建拦沙坝的作用及取水口布置等具有重要意义。

1 水库三维水沙数学模型 1.1 模型简介 1.1.1 基本方程

本文采用水动力学模型的控制方程为Reynolds时均NS方程^[15-16]，满足静压假定和Boussinesq涡黏性假定，在Cartesian坐标下，其方程表示如下：

$ \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0 $

(1)

$ \begin{aligned} \frac{\mathrm{D} u}{\mathrm{D} t}=& f v-g \frac{\partial \eta}{\partial x}-\frac{1}{\rho_{0}} \frac{\partial p_{A}}{\partial x}-\frac{g}{\rho_{0}} \int_{z}^{\eta} \frac{\partial \rho}{\partial x} d z+\\ & \frac{\partial}{\partial z}\left(K_{\mathrm{mv}} \frac{\partial u}{\partial z}\right)+K_{\mathrm{mh}}\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\right), \end{aligned} $

(2)

$ \begin{aligned} \frac{\mathrm{D} v}{\mathrm{D} t}=&-f u-g \frac{\partial \eta}{\partial y}-\frac{1}{\rho_{0}} \frac{\partial p_{\mathrm{A}}}{\partial y}-\frac{g}{\rho_{0}} \int_{z}^{\eta} \frac{\partial \rho}{\partial y} d z+\\ & \frac{\partial}{\partial z}\left(K_{\mathrm{mv}} \frac{\partial v}{\partial z}\right)+K_{\mathrm{mh}}\left(\frac{\partial^{2} v}{\partial x^{2}}+\frac{\partial^{2} v}{\partial y^{2}}\right). \end{aligned} $

(3)

其中：(x，y)表示水平Cartesian坐标，m；z表示垂向坐标，向上为正，m；(u，v，w)表示3个方向的流速，m/s；t表示时间，s；f表示Coriolis力系数，s^-1；η为自由水面水位，m；ρ₀、ρ分别表示参考密度和浑水的密度，kg/m³；g为重力加速度，m/s²；K_mh、K_mv分别表示水平与垂直涡黏性系数，m²/s；p_A表示自由水面的大气压强，N/m²。

对于自由水面采用水位函数法处理，对连续方程(1)沿水深方向积分，可得自由水面方程为

$ \frac{\partial \eta}{\partial t}+\frac{\partial}{\partial x} \int_{z_{{\rm b}}}^{\eta} u \mathrm{d} z+\frac{\partial}{\partial y} \int_{z_{{\rm b}}}^{\eta} v \mathrm{d} z=0. $

(4)

其中：z_b表示河床底高程，m。模型在河床底面的动力学边界条件由床底摩擦剪应力和水体底层的Reynolds应力平衡给出：

$ K_{\mathrm{mv}}\left(\frac{\partial u}{\partial z}, \frac{\partial v}{\partial z}\right)=\left(\tau_{b x}, \tau_{b y}\right), z=z_{\mathrm{b}}. $

(5)

其中：τ_bx、τ_by为床面的摩擦剪应力，m²/s²。在水面，不考虑风应力的影响。在开边界采用Dirichlet边界条件或Neumann边界条件，在固壁边界采用有滑移无穿透的边界条件。

1.1.2 泥沙模型

在德泽水库中的泥沙以悬移质运动形式为主，因此本研究的泥沙控制方程采用三维悬沙输运方程。三维悬沙输运方程也是基于带有沉速的对流扩散方程^[21]，

$ \frac{\partial C}{\partial t}+\frac{\partial(u C)}{\partial x}+\frac{\partial(v C)}{\partial y}+\frac{\partial\left[\left(w-\omega_{s}\right) C\right]}{\partial z}=\\ \frac{\partial}{\partial z}\left(K_{\mathrm{sv}} \frac{\partial C}{\partial z}\right)+K_{\mathrm{sh}}\left(\frac{\partial^{2} C}{\partial x^{2}}+\frac{\partial^{2} C}{\partial y^{2}}\right). $

(6)

其中：C为泥沙浓度，kg/m³；ω_s为泥沙颗粒沉速，m/s；K_sh、K_sv分别为水平和垂直方向的泥沙扩散系数，m²/s，水平扩散系数选取常数处理，垂直扩散系数根据紊流模式封闭。水面处的边界条件为

$ \left[\omega_{\mathrm{s}} C+K_{\mathrm{sv}} \frac{\partial C}{\partial z}\right]=0. $

(7)

在床面，采用通量型边界条件，具体形式为

$ \left[\omega_{{\rm s}} C+K_{\mathrm{sv}} \frac{\partial C}{\partial z}\right]=\omega_{\mathrm{s}}\left(C_{\mathrm{b}}-\alpha_{*} S_{*}\right). $

(8)

其中C_b为床面含沙量，kg/m³。可根据文[22]的方法计算

$ C_{\mathrm{b}}=C_{\mathrm{m}}+\left[1-\exp \left(-\frac{\omega_{\mathrm{s}}}{K_{\mathrm{sv}}}\left(z_{\mathrm{m}}-z_{0}\right)\right)\right] \alpha_{*} S_{*}. $

(9)

其中：C_m和z_m分别为底层网格中心处的含沙浓度和距河床底高度；α_*为非平衡输沙响应系数，表示床面含沙量与垂向平均含沙量之比^[23]，可随水沙条件动态调整；S_*为垂向平均挟沙力，kg/m³。选择在黄河上应用广泛的张红武等^[24]的公式计算：

$ S_{*}=2.5 \times\left[\frac{\left(0.0022+C_{\mathrm{V}}\right) U^{3}}{\kappa_{\mathrm{m}} \frac{\gamma_{\mathrm{s}}-\gamma_{\mathrm{m}}}{\gamma_{\mathrm{m}}} g h \omega_{\mathrm{m}}} \ln \left(\frac{h}{6 D_{50}}\right)\right]^{0.62}. $

(10)

其中：U为垂向平均流速；C_V为体积含沙量；h为水深；κ_m、ω_m、γ_s、γ_m分别为浑水Karman系数、非均匀沙群体沉速、泥沙颗粒容重及浑水容重。

与式(8)相对应的河床变形方程为

$ \rho^{\prime} \frac{\partial Z}{\partial t}=\omega_{{\rm s}}\left(C_{{\rm b}}-\alpha_{*} S_{*}\right). $

(11)

其中：ρ′为河床组成物质的干密度，kg/m³。在进口边界，采用Dirichlet边界给定悬沙的浓度和级配条件；在出口处，采用Neumann边界条件，设定泥沙浓度梯度为0。对于非均匀沙，采用李义天^[25]方法求得不同粒径组的分组挟沙力，对不同粒径分组分别进行上述过程的计算。

1.2 研究对象 1.2.1 计算区域与网格划分

德泽水库坝高142 m，正常蓄水位1 790 m，正常库容41 597万m³，最低运行水位1 752 m，相应死库容18 902万m³。模拟范围包括整个德泽水库库区，由干流牛栏江与支流干河组成，位于库区的牛栏江干流长约40 km，支流干河长约5 km，干河约在距坝址17.7 km处汇入牛栏江干流，库区平面如图 1所示。泵站取水口布置于干、支流交汇口下游约330 m处左岸，取水口断面天然河床高程1 730 m，取水口底板高程1 744 m，设计取水流量23 m³/s。

德泽水库为峡谷型水库，河道狭窄，则河道范围内采用四边形网格，网格尺度约2~10 m，其余地区采用三角形网格，网格尺度为10~40 m，共得到100 062个网格。计算初始地形条件由2017年实测地形插值得到，插值后的地形如图 1所示，干、支流交汇区域及取水口附近的地形和网格划分如图 2所示。此外，为了更好地模拟河道至坝前水深变化幅度大的特点并且提高计算效率，垂向上采用分层SZ坐标^[26]，共分45层，其中Z坐标30层，每层厚度约3 m，S坐标15层，每层厚度约0.2~1.5 m。

根据现场考察发现，干、支流均有修建拦沙坝的合理位置^[1]，如图 2所示。设计干河拦沙坝位于干河取水口断面上游2.3 km，溢流坝堰顶高程1 785 m，当水库运用水位高于1 785 m时，干河拦沙坝处于全淹没状态。孔口总宽度21 m，底孔不参与泄洪，底板高程1 760 m。设计牛栏江拦沙坝位于牛栏江取水口断面上游3.4 km处，溢流坝堰顶高程1 795 m，该拦沙坝始终不会被完全淹没。孔口总宽度42 m，底孔不参与泄洪，底板高程1 765 m。

现有取水口的高程如图 3a所示，在低水位时距回水末端较近，且垂向上靠近河底。而新建分层取水方案充分利用了取水口下游侧李子箐冲沟的有利地形新建分层取水口，通过短隧洞与现引水隧洞相接，如图 3b所示，利用围堰保证水库的表层水进入进水前池，则可从新建的取水口取到表层水。

1.2.2 边界条件与相关参数

德泽水库的进口包括牛栏江、干河入口以及区间入汇，其中区间汇入水沙的过程主要集中在德泽大坝上游约5 km处的暗河处，对取水口没有影响，本文不予讨论。德泽水库的出口包括发电引水口、泄水建筑物和泵站取水口。进口边界条件为牛栏江、干河和区间入汇的流量和含沙量过程，入库水温分别为20、17和18 ℃；出口边界条件为根据水量平衡和坝前控制水位得到的发电引水口的流量过程。泵站取水口则根据实际运行情况给定取水流量过程。库内初始水温采用20 ℃，初始悬沙含量为0。

根据现场观测资料将入库悬移质和床沙质的级配分为5组，选用的代表粒径分别为0.001、0.006、0.012、0.02和0.06 mm，各组泥沙所占的比例如表 1所示，相应干、支流中值粒径分别为0.007 2和0.005 3 mm。可根据文[27]确定各组泥沙的沉速，并考虑泥沙浓度对沉速的影响^[9]得到各代表粒径的沉速如表 1所示。

1.3 模型验证

选用2017年7月19—21、26日的牛栏江及干河现场观测结果作为第3—5、10天的入库水沙资料，插补延长得到12天验证水沙系列，计算时间步长为10 s。牛栏江和干河的水沙过程如图 4所示，并根据多年实测资料假定坝前暗河的流量为牛栏江流量的0.2倍，含沙量与牛栏江一致。采用无拦沙坝的现状地形，并根据实测情况设置运行水位在1 785 m左右。

由于实测牛栏江入库水沙资料位于模型设置入口下游约10 km处，因此干流的计算结果与实测结果没有可比性，故不对牛栏江的验证结果进行分析。选择模型计算干河第5天和第10天的沿程表、底层含沙量与实测值进行对比，如图 5所示。由图可知，模拟的干河含沙量分布与实测结果有较为一致的变化趋势，且图 5a和5b表明第5天时沿程的底层含沙量缓慢降低，而表层含沙量迅速降低，与实测的异重流现象较为一致。但采用ADCP走航测量沿程含沙量时，测量船受到水力和风力的影响，发生位置偏移，导致含沙量出现波动，而计算值无法反应此类不可控因素造成的含沙量波动，因此两者之间存在一定的误差(见图 5a)。但从两者的变化趋势和含沙量大小来看，模型计算的结果是可靠的，说明该模型能够用于模拟德泽水库的水沙运动过程。

2 方案模拟结果分析

在验证计算的基础上，为分析修建拦沙坝前、后沿程及取水口分层含沙量的变化情况，分别以无拦沙坝和有拦沙坝作为初始地形条件，选择1997年作为大水大沙年(径流量22.5亿m³、输沙量186.4万t)、1998年作为中水中沙年(径流量16.8亿m³、输沙量143.5万t)、2006年作为小水小沙年(径流量9.5亿m³、输沙量68.7万t)，以这3种典型年的汛期水沙过程为三维数模计算的水沙条件。由于11月流量及含沙量较小，本文计算的汛期过程指6—10月共153天的过程，坝前水位过程根据水库调度运用方式得到。各典型年的泵站取水流量过程一致，7月份时停泵检修，其余时间维持23 m³/s。

2.1 干、支流拦沙坝的拦沙效果

在泵站取水口以上河段修建拦沙坝后，通过壅高拦沙坝上游河段的水库运用水位、降低水流流速，从而达到促进入库泥沙在拦沙坝以上库区淤积的目的。根据式(12)可分别计算有、无拦沙坝地形条件下干、支流入库—拦沙坝河段的拦沙率(即淤积比)，从而分析拦沙坝的拦沙效果。

$ \phi_{\mathrm{aver}}=\frac{W_{\mathrm{sin}}-W_{\mathrm{sout}}}{W_{\mathrm{sin}}} \times 100 \%. $

(12)

其中：ϕ_aver表示平均拦沙率；W_sin表示干、支流的入库沙量，万t；W_sout表示拦沙坝断面的输沙量，万t。不同典型年汛期过程的拦沙率如表 2所示，由表 2可知，无拦沙坝时，干、支流入库—拦沙坝段的拦沙率基本呈现来水来沙量越大，拦沙率越小的变化特点。有拦沙坝时，干河的拦沙率在大、中、小沙年的拦沙率分别为41.46%、45.40%和49.93%，较无坝时绝对增加6.21%、2.67%和4.27%；牛栏江入库—拦沙坝的拦沙率在大、中、小沙年的拦沙率分别为51.34%、53.95%和58.02%，较无坝时绝对增加17.32%、15.05%和15.85%。可见，拦沙坝的存在能增加入库—拦沙坝库区的拦沙率，促进泥沙在拦沙坝上游淤积，减小泥沙对取水口断面的影响。

2.2 水库沿程含沙量分布特点

拦沙坝的存在对干、支流的水流运动和泥沙输移有明显影响。以中水中沙年为例，分别选取有、无拦沙坝情况下较为典型的干、支流沿程含沙量分布，如图 6-9所示。

2.2.1 干河沿程含沙量分布

图 6和7分别表示无坝和有坝条件下干河的沿程含沙量分布。由图可知，拦沙坝修建前后，泥沙均在输移过程中不断沉降，形成表层小底层大的含沙量垂线分布，符合泥沙运动的一般规律。图 7表示拦沙坝壅高上游水位，导致了水流流速降低，相同时间内洪水传播距离变短，异重流的潜入点位置向上游移动，如图 6显示第44天在距汇口约3.5 km处浑水潜入形成异重流，表层含沙量迅速降低，而图 7在第44天时异重流潜入位置距汇口约4.2 km处。对比图 6和7可知，在同一时间，有拦沙坝的坝下游河段的含沙量明显低于无坝时相同河段的含沙量，说明拦沙坝拦截了一部分泥沙，减少了进入下游的泥沙。

2.2.2 牛栏江沿程含沙量分布

图 8和9分别表示无坝和有坝条件下牛栏江的沿程含沙量分布。由图可知，拦沙坝修建前后，牛栏江沿程含沙量分布的变化情况与干河的情况基本一致。但由于牛栏江入库流量较大，含沙量较小，在牛栏江上形成异重流的机会少于干河中形成异重流的机会。此外，图 8显示在第44天干流洪水未传递到交汇口时，交汇口附近河段已经出现较高的含沙量，由于干、支流的入口至交汇口距离不同，洪水抵达交汇口的时间不同，导致交汇区域出现干、支流相互入侵的现象，对附近河段的含沙量垂向分布扰动较大，从而影响取水口的含沙量分布。

2.3 取水口分层含沙量变化特点

根据取水口断面各层的含沙量和流速可得到取水口断面的平均含沙量，不同典型年取水口断面平均含沙量的汛期平均值如表 3所示。由表可知，无拦沙坝时，大、中、小沙年的取水口断面平均含沙量分别为0.28、0.21和0.24 kg/m³，有拦沙坝时降低为0.19、0.14和0.15 kg/m³，分别减少32.9%、31.5%和36.7%，说明拦沙坝能够降低取水口断面的平均含沙量。

为了研究取水口断面分层含沙量的变化情况，以水面层、取水口高度层和河床附近层为表层、中层、底层，拦沙坝修建前后，中水中沙年取水口断面各层含沙量的变化如图 10所示。由图可知，无拦沙坝时，取水口上游河段发生了异重流现象，取水口断面含沙量分层明显，呈现表层、中层、底层含沙量由小到大的分布特点。有拦沙坝以后，拦沙坝破坏了明渠流泥沙垂线分布和异重流运动的连续性，水流过坝后泥沙在垂向上重新分布，导致取水口断面的含沙量趋于均匀分布。

为分析不同典型年取水口中、表层含沙量在有、无拦沙坝条件下随时间的变化情况，以中水中沙年条件为例，变化情况如图 11所示。由图可知，有坝时，中层含沙量在大部分时间小于无坝时的中层含沙量，而表层含沙量却大于无坝时的表层含沙量，说明拦沙坝能够减小取水口断面的中层含沙量，但是却会造成表层含沙量的增加。不同典型年的中、表层含沙量的汛期平均值如表 3所示，无坝时大、中、小沙年的中层含沙量分别为0.27、0.21和0.24 kg/m³，表层含沙量分别为0.07、0.06和0.07 kg/m³，有坝时大、中、小沙年的中层含沙量分别降低了45.7%、42.2%和41.4%，而表层含沙量却分别增加了63.5%、59.2%和34.9%。

通过对瀑布公园及盘龙江水质进行的现场调研发现^[2]，瀑布景观具有良好视觉状态可以接受的浊度上限值为17.1 NTU。现场观测含沙量与浊度TU的对应关系为

$ S=0.001156 \times \mathrm{TU}+0.03652. $

(13)

其中：S为含沙量，kg/m³，浊度的单位为NTU。由式(13)可得17.1 NTU对应的含沙量为0.056 kg/m³，将不同典型年的取水口中、表层含沙量按式(13)换算成相应的浊度值，并将取水口浊度在3个区间(< 17.1 NTU、17.1~100 NTU、>100 NTU)的时间分布情况列于表 4。

由表 4可知，无拦沙坝时的中层浊度 < 17.1 NTU的天数占汛期的比例为10.8%~16.6%，>100 NTU区间的比例为32.8%~61.6%；表层浊度 < 17.1 NTU的天数占汛期的比例为47.0%~62.1%，>100 NTU区间的比例为0~7.6%。这说明无拦沙坝时，在 < 17.1 NTU区间内的表层天数明显多于中层天数，而在>100 NTU区间内中层天数则明显多于表层，采用分层取水能显著增加取水浊度满足要求的时间。

有拦沙坝时，中层浊度 < 17.1 NTU的天数占汛期的比例为13.9%~18.5%，>100 NTU区间的比例为28.7%~38.9%，可知 < 17.1 NTU的比例较无坝时略有增加，而>100 NTU的比例明显减小，说明拦沙坝能够有效降低中层>100 NTU的浊度，但较难将其降低至 < 17.1 NTU的范围内，对提高中层水体达标时间的比例作用不大。有拦沙坝时的表层浊度 < 17.1 NTU的天数占汛期的比例为13.6%~23.0%，>100 NTU区间的比例为16.3%~26.7%，可知 < 17.1 NTU的比例较无坝时明显减小，而>100 NTU的比例明显增加，说明拦沙坝的存在将会增加表层水体的浊度，降低表层水体达标时间的比例。

3 结论

本文利用三维水沙数学模型分别模拟了不同典型年的汛期水沙过程在德泽水库中的运动规律，对比分析了修建拦沙坝前、后沿程及取水口分层含沙量的变化情况。结果表明：

1) 现状条件下取水口断面的汛期含沙量具有明显的分层特点，表层水体的浊度有47.0%~62.1%的时间达到小于17.1 NTU的要求，其余时间的浊度虽然不能达标，但仅有0%~7.6%的时间大于100 NTU。而中层水体的浊度只有10.8%~16.6%的时间小于17.1 NTU，有32.8%~61.1%的时间大于100 NTU。可见与现有中层取水口相比，新建分层取水口方案充分利用了取水口含沙量的分层特点，在提高汛期取水浊度满足要求的时间上有明显优势，可作为牛栏江滇池补水工程汛期取水防沙措施的首选方案。

2) 修建拦沙坝可促进泥沙在拦沙坝以上库区的淤积，明显降低取水口断面的汛期平均含沙量和中层含沙量，将平均含沙量从0.21~0.28 kg/m³降至0.14~0.19 kg/m³，中层水体浊度大于100 NTU的时间比例从32.8%~61.6%降低至28.7%~38.9%。但降低后的中层含沙量仍偏高，只能将小于17.1 NTU的占比从10.8%~16.6%提升到13.9%~18.5%，效果并不明显。本研究结果说明拦沙坝方案能缓解汛期水体浑浊对滇池水环境及瀑布公园和盘龙江水景观的影响，但对提高取水浊度满足要求时间的作用有限。

3) 由于拦沙坝对泥沙垂向分层的扰动，导致取水口含沙量的垂向分布趋于均匀，表层水体浊度增加，其小于17.1 NTU的时间比例从47.0%~62.1%降低至13.6%~23.0%，拦沙坝与分层取水方式之间存在矛盾，需要注意。