基于高浊度事件调查的湖泊粘性泥沙输移模型建模方法

技术领域

本发明属于水环境及数值模拟技术领域，具体涉及基于高浊度事件调查的湖泊粘性泥沙输移模型建模方法。

背景技术

湖泊表层底泥是湖泊生态系统的重要组成部分，与波浪作用有关的高浊度事件是浅水湖泊的一个特征(Chung et al.,2009；Valipour et al.,2017)。悬浮沉积物在营养盐和污染物的生物地球化学循环中扮演重要角色(Bao et al.,2020；Carpenter,2005；Scheuet al.,2015；Sterner et al., 2017)。在湖泊高浊度事件中发生的沉积物迁移过程可以概括为：当湖床受到波浪或水流引起的剪切应力时，泥沙可能被夹带(即沉积物侵蚀)；然后悬浮沉积物通过水动力过程(平流和扩散)进行迁移；随着风的消退，悬浮沉积物又沉降到沉积床(沉积)。近年来，泥沙再悬浮的过程受到了极大的关注(Gibbs et al.,2016；Schwabet al.,2006；Troud et al.,2019)。然而，与已经研究了两个多世纪的河道环境(河流)中的沉积物迁移相比(

湖床沉积物是水、粘土、有机物、淤泥和沙子的混合物，当细粒成分质量(粒径小于62.5μm) 的百分比大于10％时，通常被称为粘性泥沙，这是沉积物运输行为从无粘性到粘性变化的临界值(Mitchener and Torfs,1996)。细粒度的成分可以大大增加床面泥沙的抗剪应力，并影响床面沉积物侵蚀的开端0(Mitchener and Torfs,1996；

由于粘性泥沙运动是高度依赖于环境条件的，因此原位研究是弄清湖泊粘性泥沙运动机理的根本方法。已有的湖泊粘性泥沙动力学研究以实验室研究为主，但是基于实验室的方法 invariably给出不准确的现场估计(Douglas et al.,2010；Liu et al.,2019；Mitchener and Torfs, 1996)。一方面，在实验室中实现与自然水体一致的水动力重建、絮凝物形成和化学属性是非常困难的；另外一方面，从采集、运输和测量过程将无可避免对粘性泥沙的质地、化学和生物条件产生影响。因此，许多泥沙学者均强烈推荐采用原位观测方法研究粘性泥沙运动 (Douglas et al.,2010；Li and Amos,2001；Scheu et al.,2015；Valipour et al.,2017)。不过，由于风场是决定湖泊粘性泥沙运动的关键因素，湖泊中显著的高浊度事件多与强风事件相关，而两者均具有间歇性和快速变化的性质。这对原位观测湖泊粘性泥沙动态提出了挑战。幸运的是，近年来的无线传感器和通信技术的发展为大风事件中湖泊环境要素变化原位观测铺平了道路 (Hipsey et al.,2019；Soulignac etal.,2017)。

对机理认识不清也导致了湖泊粘性泥沙输移模型研究的滞后。相对于具有明确机制的非粘性泥沙输移模型，已有的粘性泥沙输移模型对描述泥沙运动的关键过程的表达往往是理想化的和经验性的。例如，现有的粘性泥沙模型将泥沙垂向运移分成沉降、稳定和侵蚀三个独立的状态，三个状态发生时段互不重叠(Chao,et al.,2008；Li and Amos,2001；Morales-Marin,et al.,2017)。但是，受自然风场不稳定性和自然泥沙粒径组成的复杂性影响，自然水体中泥沙运动的三种状态并非是可以区分的独立状态。此外，由于现场数据的不足，大部分湖泊粘性泥沙输移模型均采用率定或者引用文献值的方式确定关键过程或参数(Cardenas et al.,2005； Chao,et al.,2008；James,et al.,1997；Liu,et al.,2019；Wang et al.,2013)，忽视了这些过程或参数高度的时空变化特性。这样做显然不能构建高精度的数值模型，也无助于加深我们对粘性泥沙运动特征的理解。

缺乏高时空分辨率的现场数据意味着对粘性泥沙动力学的了解仍然很有限，阻碍了大型浅水湖泊粘性泥沙模型的发展。

发明内容

为克服上述现有技术存在的问题，本发明提供一种基于高浊度事件调查的湖泊粘性泥沙输移模型建模方法。

为实现上述技术目的，本发明采用如下技术方案：

基于高浊度事件调查的湖泊粘性泥沙输移模型建模方法，包括：

(1)在湖中设置多个沉积物采样点、原位观测点；

(2)在开始原位观测之前，在原位观测点处采集沉积物柱，测定沉积物理化特性，包括中值粒径、新鲜沉积物密度、含水率和悬浮物沉降速率；

(3)收集开始原位观测之前某一年度所有沉积物采样点的表层沉积物中值粒径和新鲜沉积物密度调查数据集，用于计算湖泊临界侵蚀切应力τ

(4)开展原位观测，观测记录原位观测点处气象、流速、波浪、浊度、悬浮泥沙浓度数据；

(5)建立湖泊三维波流耦合数值模型WCCM，获取x、y、σ方向的三维水流速度u、v、w，水深H、水密度ρ、波高H

(6)建立湖泊三维波流耦合数值模型WCCM和泥沙传输模型耦合的三维粘性泥沙输运模型WCCM-CSTM，其控制方程如下：

式中，S为悬浮泥沙浓度；w

(7)基于(2)～(6)获取的数据求解WCCM-CSTM，其中，w

作为一种优选的实施方式，所述(1)中，采集的沉积物柱取表层1cm的样品中的部分测定中值粒径、新鲜沉积物密度、含水率，其余部分用于泥沙静态沉降实验，以确定悬浮物沉降速率。

作为一种优选的实施方式，所述湖床悬浮泥沙通量J

J

作为一种优选的实施方式，泥沙沉积通量D和泥沙再悬浮通量E的数值基于悬浮泥沙沉降速率w

作为一种优选的实施方式，控制方程求解时，对流相关项采用迎风差分格式求解，扩散相关项采用中心差分格式求解，垂直扩散项采用交替方向隐式差分格式求解。

作为一种优选的实施方式，所述方程求解时，使用步骤(3)获取的数据集计算湖泊临界侵蚀切应力τ

假设τ

对原位观测点的τ

将步骤(3)获取的数据集代入τ

作为一种优选的实施方式，所述(7)中，w

作为一种优选的实施方式，模型求解时，使用低参考侵蚀率值和高泥沙侵蚀系数值以描述高浊度事件中悬浮泥沙浓度的急剧增加；其中，参考侵蚀率值和泥沙侵蚀系数值通过(4) 中原位观测过程获取。

作为一种优选的实施方式，还包括，基于实测数据和/或遥感反演数据对模型计算数据进行验证。

本发明在高浊度事件综合调查的基础上，分析获取大型浅水湖的粘性泥沙输移过程情况，并根据沉积物调查、实验室实验和现场观测，得出粘性泥沙再悬浮、运移和沉积过程的数学表达式，利用这些数学表达式，在波浪-海流耦合模型的基础上建立了三维粘性泥沙输运模型并通过与来自现场观测和卫星图像检索的高频浊度的比较，评估了WCCM-CSTM的性能，验证了其应用，本发明的方法结合系统性现场认知，改进了对浅水湖泊中风生水动力驱动的高浊度事件的认识，建立了更符合实际的粘性泥沙悬浮和沉降过程数学方程，并强调了系统性现场研究对研发高精度的湖泊粘性泥沙输移模型的必要性。

附图说明

图1是太湖湖底高程等值线和东山气象站、原位观测点(F1-F5)、沉积物采样点和水位站(WL1-WL5)示意图。

图2是F1、F2站床层沉积物粒径。

图3是泥沙沉降实验中测量和计算的泥沙沉降速率。

图4是2018年野外观测风速风向变化。

图5是2018年原位观测期间F1站(A图)和F2站(B图)处实测湖流和风浪引起的剪应力变化。

图6是2018年原位观测期间F1站(A图)和F2站(B图)处实测浊度和悬浮泥沙浓度(S)变化。

图7是2018年原位观测期间据方程(19)计算得到的F1站(A图)和F2站(B图)处悬浮泥沙浓度(S)与由实测浊度转换得到的F1站(A图)和F2站(B图)处悬浮泥沙浓度对比。

图8是2020年原位观测期间风速和风向变化过程。

图9是2020年原位观测期间WCCM-CSTM计算得到的和实测浊度转换得到的F2(A图)， F3(B图)，F4(C图)和F5站(D图)处悬浮物浓度对比。

图10是2020年原位观测期间WCCM-CSTM计算得到的F2(A图)，F3(B图)，F4(C 图)和F5站(D图)处悬浮泥沙浓度垂直剖面时间变化过程。

图11是2020年12月16日由MODIS卫星影像提取的悬浮泥沙浓度场与由WCCM-CSTM计算的悬浮泥沙浓度场对比。

具体实施方式

下面以太湖为例，对本申请基于高浊度事件调查的湖泊粘性泥沙输移模型建模方法做具体阐述。

实施例中，中值粒径d

为了评估WCCM-CSTM的性能，考虑了p<0.05和p<0.01两个显著性水平下测量值和模拟值之间的相关系数(r)。使用SPSS18.0软件(SPSS Inc.,Chicago,IL,USA)计算回归分析和双变量相关性。此外，使用ArcGIS 10.2(美国ESRI Inc.)处理空间数据，使用Tecplot360 (美国Tecplot Inc.)绘制SSC等值线。

实施例1

位于长江三角洲平原区的太湖多年平均水深1.9m，水面面积2338.1km

(1)在湖中设置116个沉积物采样点、5个原位观测点(F1-F5)，如图1所示，在太湖南部湖区(F1)和湖心区(F2)设置2个原位观测点；在全湖设置4个常规原位观测点(F3-F5)；

(2)在2018年原位观测开始之前，考察F1和F2点表层底泥理化特性，采用柱状采泥器在F1和F2点处各采集3根沉积物柱。取6根沉积物柱的表层1cm样品的部分用于测定中值粒径(d

F1站和F2站沉积物d

泥沙静态沉降实验：

取F1和F2站点处表层1cm底泥风干样品配制450mg/L的悬浊液。将搅拌均匀的悬浊液注入内径10cm，高度54cm，底面积314.16cm

泥沙静态沉降实验表明(图3)：F1站沉积物实测w

此外，收集2014年全湖116个沉积物采样点(图1)的表层1cm沉积物中值粒径d

(3)为了减少漂浮蓝藻颗粒对观测结果的影响，选择2018年和2020年冬季对太湖沉积物悬浮和沉降过程开展两次原位观测。2018年12月24日18：00至2019年1月1日18：00间，对F1、F2两站处开展现场观测，记录同步观测的气温、风速、风向、相对湿度、气压、云量、水位、水深(H,m)、流速、波高(H

2020年11月23日0：00至2021年1月16日0：00间，同步记录F2-F5原位观测站处浊度变化过程(图1)。同时记录气温、风速、风向、相对湿度、气压、云量和水位，采集的数据统称为2020年现场观测数据，用于评估WCCM-CSTM的性能。

位于东山半岛的东山气象站提供了2018年和2020年原位观测的每小时气温、风速、风向、相对湿度和云量。

2018年实地观测时，风速范围为0.1～13.9m s

2020年野外观测时，风速范围为0～19m s

在湖床上安装一个仰视的声学多普勒电流剖面仪(XR；Sontek Inc.，美国)用于记录湖流，每30分钟以1Hz的频率采集以获得沿水柱的1200s整体平均值。空白和单元格大小分别为0.5m和0.4m。

将MIDAS DWR定向波形记录仪(Valeport Ltd.，UK)安装在湖床上，通过压力传感器记录PUV数据，并使用电磁电流传感器记录流速。MIDAS DWR记录的压力和流速分辨率分别为0.005m和0.001m s

YSI 6600探头带有擦拭式浊度传感器(Yellow Springs Instruments，美国)用于记录水面以下0.5m处的浊度变化，采样间隔为30分钟。YSI 6600探测器每周校准一次。此外，采用可通过客户端软件远程控制的自动水样采样系统(FC-9624；中国北京抓手科技有限公司)采集中间水层的水样，采样间隔为30分钟。采集到的水样暂时存放在充满液氮的FC-9624样品储存箱中(储存时间<6h)。水样用于测定悬浮泥沙浓度。

此外，WL1-WL5记录的水位由中国水利部提供。

(4)建立三维粘性泥沙输运模型(WCCM-CSTM)：

①控制方程；

在三维波流耦合数值模型(WCCM；Wu et al.,2021)的基础上，建立三维粘性泥沙输运模型。WCCM为三维水流模型和波浪模型耦合模型，可输出x、y、σ方向的三维水流速度(u、 v、w；m s

垂向压缩坐标系(x,y,σ,t)中的描述泥沙传输方程(即WCCM-CSTM控制方程)为：

式中，S为悬浮泥沙(S)浓度(g m

其他参数如表1所示。

表1 WCCM-CSTM中悬浮泥沙输运的参数值和变量方程

垂直扩散系数D

l和R

其中κ为冯卡门常数，z

②边界条件

湖面悬浮泥沙通量为0(σ＝1)，湖床悬浮泥沙通量(J

J

其中，D由w

其中，τ为湖流切应力(τ

其中，u

从已发表的描述泥沙沉降速率公式来看，w

w

w

Partheniades(1965)提出的等式用于计算泥沙侵蚀通量：

τ

τ

③方程求解

对流相关项采用迎风差分格式求解，扩散相关项采用中心差分格式求解，垂直扩散项采用交替方向隐式差分格式求解。

用Ji(2008)提出的一种简单的方法来模拟沉积物床的变化，方程如下：

④模型设置

为消除水动力初始值和SS浓度对2019年野外观测模拟的影响，计算时间相对2020年野外观测提前一周。WCCM-CSTM的计算网格包括5184个水平单元和五个垂直层。

将2018年和2020年野外观测的气象数据输入模型。由于太湖在模拟中可以认为是一个封闭的湖泊，因此可忽略流入和流出的沉积物。湖床设置为无限的沉积物来源(Cardenas等， 2005；Lou等，2000)。初始悬浮泥沙浓度、水位、流速分别设置为20mg L

WCCM-CSTM的水动力参数基于现有文献确定(Wu等，2021)。沉积物参数(表1) 基于高浊度事件的调查确定。

已有的泥沙沉降实验表明太湖泥沙平均沉降速度不超过2.0×10

相关性分析表明(图3)，使用方程(16)计算的w

表2大型湖泊粘性泥沙输运研究中w

图4为2018年野外观测风速风向变化，图5为2018年原位观测期间F1站(A图)和 F2站(B图)处湖流和风浪引起的剪应力变化观测结果。F1和F2站的当前速度范围分别为 0.1-12.8 cm s

表3 2018年原位观测期间F1和F2站处实测风速、水流切应力、波浪切应力、浊度相关性分析

**表示在p<0.01水平上相关性显著(双尾检验).

F1站浊度从57.8(12月26日3:00)急剧上升至231.4 NTU(12月27日16:00)，12月29日3:00后逐渐下降。悬浮泥沙浓度(S)与浊度呈显著正相关(r＝0.93，p<0.01；图6A)。关系为：

S＝2.05Turbidity-53.9 (17)

同时，F2站浊度从35.7(12月26日4：30)急剧上升到157.5 NTU(12月27日17:00)，然后保持高浊度状态。F2站测得的S值范围为86.3～297.1mg L

S＝1.46Turbidity+23.2 (18)

浅水湖泊中，任何地点的悬浮泥沙浓度(S)的演变都可以由垂向泥沙通量J

S

根据Luettich等人提出的方法(1990)，删除方程(1)中的水平对流和扩散项，利用垂直沉积物通量估算F1和F2站的S。随后基于w

其中，F1和F2站的τ

方程(19)～(22)计算结果与2018年实地观测期间F1和F2站的S实测结果比较如图7所示，同时，根据校准后的τ

τ

利用2014年高空间分辨率沉积物调查数据集计算太湖τ

(5)WCCM-CSTM的验证

利用方程(1)-(6)、(20)-(22)开发的WCCM-CSTM用于模拟2020年实地观察期间通过1.2倍浊度计算的S(图9)。结果表明WCCM-CSTM模拟的SSC在F2、F3、F4 或F5站与实测值显著正相关(r≥0.9,p<0.01；图9)。WCCM-CSTM可以准确模拟每次高浊度事件期间SSC的急剧增加。此外，任意两个站点之间的SSC也显著相关(p<0.01)。每条S曲线(图9)主要由五次高浊度事件组成。各站浊度最高的事件发生分别在2020年12 月29日和2021年1月3日。

WCCM-CSTM计算的F2～F5站处悬浮泥沙浓度垂直剖面随时间变化显示了发生在2020 年12月29日和2021年1月3日间的高浊度事件的发生过程(图10)。在这四个原位观测站点中，F5站处的S浓度垂向浓度变化梯度值从高浊度事件前的0.06mg L

太湖S是从美国国家航空航天局(http://modis.gsfc.nasa.gov/)2020年11月16日至2021年 1月23日的中分辨率成像光谱仪(MODIS,250m)1b级数据中提取的。首先对MODIS数据进行大气校正算法处理，然后采用为太湖开发的S定量反演模型(Jiang et al.,2009)利用MODIS 数据反演S。

在五次高浊度事件中，利用2020年12月16日的少云或无云的MODIS图像反演S场(图 11)，并与WCCM-CSTM模拟进行比较表面S场。结果表明，MODIS反演S场的高值区域主要分布在西南湖，这与相应的WCCM-CSTM模拟S场的空间分布相似。

太湖是一个大型浅水湖泊，受风生水体运动影响，频繁发生高浊度事件。沉积物调查表明，太湖表层底泥d

风浪是引发太湖高浊度事件的主要驱动力。2018年的原位观测期间，F1和F2站的τ

高浊度事件可由侵蚀过程和沉降过程描述，这两个过程是可以同时发生的。Luettich et al. (1990)提出浅水湖泊本地的泥沙垂直通量决定高浊度事件，水平输送的贡献可以忽略。以此理论，τ

太湖粘性泥沙沉降过程中的关键参数w

描述太湖底泥侵蚀过程的关键参数为M和τ

尽管太湖浊度与S比例具有较大的时空差异，但是WCCM-CSTM能够精确模拟水体SSC 的变化趋势。2020年现场观测的54天内，WCCM-CSTM模拟的F2、F3、F5和F6站水体的 SSC与的实测SSC相关系数在0.89以上，五次浊度事件的最高值及其发生时间也较为吻合(图9)。模拟值与实测值的绝对误差源于浊度与S换算。2018年原位观测显示，F1站S与浊度的回归方程不同于F2站(图6)，也不同于其他太湖观测研究得到的回归方程(Ding et al.,2012；Tang et al.,2019；Zheng et al.,2015)。同时，每个站点S与浊度的比例系数在时间上也存在较大变化。这是因为作为一个遭受严重蓝藻水华影响的湖泊，太湖浊度不仅受底泥悬浮的影响，还受到水体中藻类颗粒物含量的影响。即使在冬季，太湖水体中仍有较高蓝藻生物量(Ma et al.,2016)。此处采用1.2的乘数是为了保证模拟值与实测值之间的绝对误差最小。

WCCM-CSTM模拟SSC场比MODIS反演S场更为可靠反映由底泥侵蚀和悬浮物沉降导致的水体SSC空间分布。Puls等(1994)认为数值模型不能成功再现通过遥感反演技术获得的泥沙空间分布。本次研究的WCCM-CSTM模拟的S场与MODIS反演SSC场存在较大的差别(图11)。主要原因为：(1)我们假定MODIS影像拍摄时间均为上午10：30，这可能与实际时间有差别；(2)太湖水体中的S随水深增大(Zheng et al.,2015)，WCCM-CSTM 模拟的S场是表层水体S均值，这显然要大于MODIS反演的水体表面的S；(3)太湖水体中的悬浮物粒径也随水深而增大(Zheng et al.,2015)，水表的悬浮物以沉降缓慢的或漂浮态的颗粒物为主，其时空运移特征与下层水体中源自底泥悬浮的颗粒物的时空运移特征有明显区别，例如漂浮蓝藻易于在下风区沿岸带聚集(Wu et al.,2019)。实际上，尽管2020年原位观测中各站间S浓度时间序列具有显著相关性(p<0.01)，但是其五次高浊度事件中的S浓度峰值存在较大差异(图9)。这种差异主要取决于各站间与底泥理化性质相关的τ

太湖高浊度事件综合调查结果表明，太湖表层沉积物以粒径小于62.5μm的粘性泥沙为主。风浪是导致高浊度事件中S浓度快速变化的主要因素。水平传输对湖心区高浊度事件的贡献要大于岸边带。高浊度事件可以分为侵蚀过程和沉降过程，但是两个过程是可以同时发生的。沉降方程中的w