1.本发明涉及水环境模拟技术领域，更为具体地，涉及一种流域水环境模 拟方法及装置。
背景技术：2.随着城市化进程和经济快速发展，我国的水污染问题比较突出，水污染 问题直接威胁人们的身体健康和生命安全。对流域水环境的综合治理得到各 级管理部门的高度重视。流域水环境模拟将为水环境治理方案设计和实施提 供科学依据。
3.作为流域水环境解析与模拟的重要手段，基于数值模型的流域水环境模 拟及其应用成为相关研究的热点。根据水环境模型的不同功能用途，可将水 环境模型分为陆地模型、水体模型以及陆地模型和水体模型的耦合模型。陆 地模型仅计算非点源污染累积和冲刷量，获得陆地污染负荷在时间和空间的 平均值(emc，event mean concentration)，但缺少陆地污染负荷随降雨汇流运 动过程，即缺少陆地污染物的输移模块，应用比较广泛的陆地模型有swmm、 swat、hspf等，但这些模型因缺少陆地污染负荷随水文汇流的输移模块， 不能预测陆地污染负荷详细的时空变化信息，仅能与水体模型进行单向耦合； 而水体模型包括水动力模型和水质模型，主要用来模拟污染物在河流、湖泊、 水库、河口、沿海等水体中的对流扩散过程，得到污染物在水体中的详细分 布信息，应用比较广泛的水体模型包括wasp、efdc、mike系列等，但这 些模型缺少陆地污染负荷进入水体的机理分析及水体过程与陆地过程的互馈 机制。基于单一的陆地模型或水体模型可以模拟流域出口的水文要素、营养 负荷及其在水体中的对流扩散等过程，割裂了流域水—陆水文、水质响应过 程。
4.与单一的陆地或水体模型相比，陆地与水体动态双向耦合模型通过耦合 陆地和水体模型，达到陆地及河湖水体污染物的综合模拟。陆地模型与水体 模型的耦合方式主要分为外部耦合和内部耦合。外部耦合中，陆地模型和水 体模型独立运行，陆地模型的输出作为水体模型的输入，驱动水体模型运行， 是单向耦合(one-way coupling)。这类耦合模型研制简便，计算速度快，适 合大尺度的流域水环境模拟。但是该方法的耦合边界点数量有限，耦合边界 点不随河湖、水库的水位涨落而变化，边界点的流量及污染物的负荷量可能 被夸大。内部耦合中，陆地模型与水体模型在时间上同步运行，陆地模型的 结果为水体模型提供边界条件，水体模型的结果通过边界条件影响陆地模型 的模拟，以此实现陆地模型和水体模型的动态双向耦合，但由于内部耦合模 式的耦合边界是固定的，当陆地模型与水体模型的网格不匹配时，固定耦合 边界点上的水位、流量、污染物变量需要通过插值实现，不能保证通过耦合 边界的通量守恒。
技术实现要素：5.鉴于上述问题，本发明的目的是提供一种流域水环境模拟方法及装置， 克服采用单一的陆地模型或水体模型模拟流域水环境污染中割裂了流域水-陆 水文水质响应过程的不足，克服目前陆地模型与水体模型耦合中，耦合边界 不明确、耦合界面通量不守恒等
问题。
6.本发明提供一种流域水环境模拟方法，包括如下步骤：
7.s1、根据预设水深阈值，对将要预测流域进行网格划分，将水深低于预 设水深阈值的网格作为非淹没网格，水深大于所述预设水深阈值的网格作为 淹没网格；
8.s2、以所述非淹没网格所形成的非淹没区与所述淹没网格所形成的淹没 区之间的边界为耦合动边界连接陆地模型和水体模型；
9.s3、对于所述非淹没网格采用所述陆地模型模拟污染物在陆地上的产生 及随径流的输移过程，计算所述非淹没区的污染物负荷分布数据；对于所述 淹没网格将所述非淹没区的污染物负荷分布数据作为计算条件，采用所述水 体模型模拟污染物在水体中的对流扩散过程，计算所述淹没区的污染物的浓 度分布数据；
10.s4、基于模型计算的时间步长，计算所述淹没网格和所述非淹没网格的 水深数据，并根据所述水深数据对所述耦合动边界进行更新，返回步骤s2。
11.此外，优选的方案是，根据预设水深阈值，对将要预测流域进行网格划 分，将水深低于预设水深阈值的网格作为非淹没网格，水深大于所述预设水 深阈值的网格作为淹没网格包括：
12.将所述将要预测流域划分为均匀的矩形网格；
13.对所述矩形网格的边界以及相邻的所述矩形网格之间的关系进行定义， 并基于所述定义对所述矩形网格进行编号；
14.将变量储存至编号的矩形网格的中心位置，给每个编号的矩形网格设置 初始条件，并将水深低于所述预设水深阈值的编号的矩形网格作为非淹没网 格，水深大于所述预设水深阈值的编号的矩形网格作为淹没网格。
15.此外，优选的方案是，所述陆地模型包括：
16.用于计算所述非淹没网格的汇流速度和水深的水文模型、用于根据所述 非淹没网格的汇流速度和水深，计算所述非淹没网格中的污染物的累积量和 冲刷量的非点源污染负荷模型和用于将所述非淹没网格中的污染物的累积量 和冲刷量作为源项，获得非点源污染物负荷分布数据的污染物输移模型；
17.所述水体模型包括：
18.用于将所述非淹没网格的汇流速度和水深作为水体水流运动的驱动条 件，通过求解局部区域的二维浅水方程，获得所述淹没网格内的流速和水深 的水动力模型和用于将所述非淹没网格的非点源污染物的浓度分布数据作为 水体中污染物对流扩散过程的驱动条件，并根据所述淹没网格内的流速和水 深，求解污染物在水体中的对流扩散过程的水质模型。
19.此外，优选的方案是，在对于所述非淹没网格采用所述陆地模型模拟污 染物在陆地上的产生及随径流的输移过程，计算所述非淹没区的污染物负荷 分布数据的过程中，
20.在所述水文模型中：
21.通过反距离插值法计算所述非淹没网格流域的降雨量，通过格林-安普顿 公式计算流域的下渗量；
22.根据所述流域的降雨量和下渗量计算所述非淹没网格的净雨量；
23.采用扩散波方程计算所述非淹没网格的汇流速度和水深；
24.在所述非点源污染负荷模型中，基于所述水文模型计算得到的汇流速度 和水深，采用指数方程计算随降雨径流进入水体的污染物的量；
25.在所述污染物输移模型中：
26.采用水质方程模拟污染物在陆地上随水流的输移过程，采用godunov中 心格式的有限体积算法对水质方程进行离散，得到所述非淹没区的污染物负 荷分布数据。
27.此外，优选的方案是，在对于所述淹没网格采用所述水体模型模拟污染 物在水体中的对流扩散过程，并将所述非淹没区的污染物负荷分布数据作为 计算条件，计算所述淹没区的污染物的浓度分布数据的过程中，
28.所述水动力模型与所述水质模型采用紧密耦合求解的方法，进行耦合求 解，采用muscl线性插值的方法将所述淹没网格中心点处的变量插值到所 述淹没网格的界面上；
29.采用godunov格式的有限体积法进行空间离散，通过hllc近似黎曼求 解器求解所述淹没网格界面上的对流通量，采用二阶runge-kutta法保证时间 积分的二阶精度，采用muscl格式保证空间上的二阶精度，求解得到所述 淹没区网格界面的水深和流速；
30.基于所述耦合动边界界面的质量、动量守恒，将所述非淹没区的污染物 负荷分布作为水体模型的计算条件，并将所述淹没网格的水深和流速作为已 知量代入到水质方程中，采用godunov格式进行空间离散，采用hllc近似 黎曼解求解通过所述淹没网格界面的污染物质的通量，采用二阶runge-kutta 法保证时间积分的二阶精度，采用muscl格式保证空间上的二阶精度，得 到所述淹没区的污染物的浓度分布数据。
31.此外，优选的方案是，所述模型计算的时间步长的确定方法包括：
32.采用显式格式的网格离散方法，确定模型计算的时间步长；其中，所述 时间步长的限制条件为：
33.其中，
34.式中，c是常数，用于保证格式的稳定性，1/2≤c＜1；u
ni,j
,v
ni,j
分别为当 前单元在x和y方向的流速分量(m
·
s-1
)；h
ni,j
为单元网格的水深(m)；g为 重力速度；d
x
和dy分别为x和y方向的扩散系数(m2·
s-1
)；γ为扩散项稳定性 系数。
35.此外，优选的方案是，基于模型计算的时间步长，计算所述淹没网格和 所述非淹没网格的水深数据，并根据所述水深数据对所述耦合动边界进行更 新，包括：
36.计算所述淹没网格和所述非淹没网格间断处的特征波，并根据所述特征 波的方向确定所述淹没区和所述非淹没区之间的质量和动量传递方向；
37.根据所述淹没区和所述非淹没区之间的质量和动量传递方向，更新所述 非淹没区和所述淹没区的水深和流速；
38.根据所述非淹没区和所述淹没区的水深，将水深低于所述预设水深阈值 的网格作为更新的非淹没网格，水深大于所述预设水深阈值的网格作为更新 的淹没网格；
39.根据所述更新的非淹没网格所形成的新的非淹没区和所述更新的淹没网 格所形成的新的淹没区之间的边界更新连接所述陆地模型与所述水体模型的 耦合动边界。
40.此外，优选的方案是，所述计算所述淹没网格和所述非淹没网格间断处 的特征波，并根据所述特征波的方向确定所述淹没区和所述非淹没区之间的 质量和动量传递方
向包括：
41.当所述耦合动边界处的特征波的方向全部指向所述淹没网格，则根据所 述非淹没网格的水深和流速计算所述淹没网格与所述非淹没网格的间断处的 通量，并判定所述非淹没区向所述淹没区传递质量，不传递动量；
42.当所述耦合动边界处的特征波在所述淹没网格和所述非淹没网格均有分 布，则根据所述淹没网格的水深及流速和所述非淹没网格的水深及流速计算 所述淹没网格与所述非淹没网格的间断处的通量，并判定所述非淹没区向所 述淹没区传递质量，所述淹没区向所述非淹没区传递动量；
43.当所述耦合动边界处的特征波的方向全部指向所述非淹没网格，则根据 所述淹没网格的水深和流速计算所述淹没网格与所述非淹没网格的间断处的 通量，并判定所述淹没区向所述非淹没区传递质量和动量。
44.本发明提供的流域水环境模拟装置，包括：
45.网格划分模块，用于根据预设水深阈值，对将要预测流域进行网格划分， 将水深低于预设水深阈值的网格作为非淹没网格，水深大于所述预设水深阈 值的网格作为淹没网格；
46.模型建立模块，用于以所述非淹没网格所形成的非淹没区与所述淹没网 格所形成的淹没区之间的边界为耦合动边界，建立陆地模型和水体模型；
47.计算模块，用于对于所述非淹没网格采用所述陆地模型模拟污染物在陆 地上的产生及随径流的输移过程，计算所述非淹没区的污染物负荷分布数据； 对于所述淹没网格将所述非淹没区的污染物负荷分布数据作为计算条件，采 用所述水体模型模拟污染物在水体中的对流扩散过程，计算所述淹没区的污 染物的浓度分布数据；
48.网格更新模块，用于基于模型计算的时间步长，计算所述淹没网格和所 述非淹没网格的水深数据，并根据所述水深数据对所述耦合动边界进行更新。
49.从上面的技术方案可知，本发明提供的流域水环境模拟方法及装置，通 过将陆地模型与水体模型直接动态双向耦合的方式，克服了采用单一的陆地 模型或水体模型以及陆地模型与水体模型单向耦合模拟流域水环境污染的不 足。通过根据预设水深阈值判断网格的状态，将水深较深的区域作为淹没区， 发生污染物的对流扩散过程，采用水体模型计算淹没区水深和流速，并求解 污染物的浓度，将水深较浅的区域作为非淹没区，发生污染物的输移过程， 采用陆地模型模拟污染物在陆地上的产生及输移过程，建立的两个模型的动 态双向耦合，计算中不指定边界位置，污染物和水流进入淹没区的地点是动 态变化的，洪水传播也可以淹没原有的汇流点，边界推移形成新的耦合边界， 实现了两种区域的交替转换，模拟更接近真实情境。对比外部的单一耦合方 式需要先后计算陆地模型和水体模型，本发明采用的双向耦合方法中两种模 型同时计算，污染物质可以由非淹没区汇入淹没区之中，而河道湖泊中的污 染物质也可以侵入汇水陆地区域，由此实现了岸上和水体中的污染物的相互 影响，更符合物理实际，更具科学性。耦合动边界作为控制体单元的边界， 在动边界两侧采用不同的计算模型(陆地模型和水体模型)，通过水量、动量、 污染物质量守恒的有限体积格式，既保证了各种物理量的守恒，也可自动实 现陆地污染物质通过耦合动边界进入水体，实现陆地模型对水体模型的驱动， 而不是通过人为地设置边界条件，将陆地污染物质作为边界条件加载到水体 模型中。
50.为了实现上述以及相关目的，本发明的一个或多个方面包括后面将详细 说明的特征。下面的说明以及附图详细说明了本发明的某些示例性方面。然 而，这些方面指示的仅仅是可使用本发明的原理的各种方式中的一些方式。 此外，本发明旨在包括所有这些方面以及它们的等同物。
附图说明
51.通过参考以下结合附图的说明，并且随着对本发明的更全面理解，本发 明的其它目的及结果将更加明白及易于理解。在附图中：
52.图1为根据本发明实施例的流域水环境模拟方法的流程示意图；
53.图2为根据本发明实施例的耦合动边界界面的移动过程示意图；
54.图3为根据本发明实施例的流域水环境模型计算流程图；
55.图4为根据本发明实施例的网格边界及相邻网格间的关系定义示意图；
56.图5为根据本发明实施例的陆地污染物负荷计算流程示意图；
57.图6为根据本发明实施例的muscl线性插值格式示意图；
58.图7为根据本发明实施例的x方向网格界面左右两侧特征线示意图；
59.图8为根据本发明实施例的y方向界面左右两侧特征线示意图；
60.图9为根据本发明实施例的二维坡面流域地形分布图；
61.图10为根据本发明实施例的二维坡面流域网格离散图；
62.图11为根据本发明实施例的初始时刻污染物分布示意图；
63.图12为根据本发明实施例的不同模型模拟的流域出口的流量过程示意 图；
64.图13为根据本发明实施例的不同时刻污染物的浓度分布示意图 (a.t＝4200s；b.t＝6300s)；
65.图14为根据本发明实施例的污染物质量守恒统计示意图。
66.在附图中相同的标号指示相似或相应的特征或功能。
具体实施方式
67.在下面的描述中，出于说明的目的，为了提供对一个或多个实施例的全 面理解，阐述了许多具体细节。然而，很明显，也可以在没有这些具体细节 的情况下实现这些实施例。
68.针对前述提出的目前对水环境模拟技术采用的陆地模型与水体模型单向 耦合的方法，忽略了水体过程对陆地过程产生的影响，割裂了流域水—陆水 文、水质响应过程，不够科学，计算精度低等问题，提出了一种动态双向耦 合的流域水环境模拟方法及装置。
69.以下将结合附图对本发明的具体实施例进行详细描述。
70.为了说明本发明提供的流域水环境模拟方法及装置，图1示出了根据本 发明实施例的流域水环境模拟方法的流程；图2示出了根据本发明实施例的 耦合动边界界面的移动过程；图3示出了根据本发明实施例的本发明实施例 的流域水环境模型计算流程；图4示出了根据本发明实施例的网格边界及相 邻网格间的关系定义；图5示出了根据本发明实施例的陆地污染物负荷计算 流程；图6示出了根据本发明实施例的muscl线性插值格式；图7示出了 根据本发明实施例的x方向网格界面左右两侧特征线；图8示出了根据本发 明实施例
的y方向界面左右两侧特征线；图9示出了根据本发明实施例的二 维坡面流域地形分布；图10示出了根据本发明实施例的二维坡面流域网格离 散；图11示出了根据本发明实施例的初始时刻污染物分布；图12示出了根 据本发明实施例的不同模型模拟的流域出口的流量过程；图13示出了根据本 发明实施例的不同时刻污染物的浓度分布(a.t＝4200s；b.t＝6300s)；图14示出 了根据本发明实施例的污染物质量守恒统计。
71.如图1所示，本发明提供的流域水环境模拟方法，包括如下步骤：
72.s1、根据预设水深阈值，对将要预测流域进行网格划分，将水深低于预 设水深阈值的网格作为非淹没网格，水深大于预设水深阈值的网格作为淹没 网格；
73.s2、以非淹没网格所形成的非淹没区与淹没网格所形成的淹没区之间的 边界为耦合动边界连接陆地模型和水体模型；
74.s3、对于非淹没网格采用陆地模型模拟污染物在陆地上的产生及随径流 的输移过程，计算非淹没区的污染物负荷分布数据；对于淹没网格将非淹没 区的污染物负荷分布数据作为计算条件，采用水体模型模拟污染物在水体中 的对流扩散过程，计算淹没区的污染物的浓度分布数据；
75.s4、基于模型计算的时间步长，计算淹没网格和非淹没网格的水深数据， 并根据水深数据对耦合动边界进行更新，返回步骤s2。
76.其中，在步骤s2建立陆地模型和水体模型的过程中，优选将陆地模型建 立在非淹没区，水体模型建立在淹没区中，以便于后续采用陆地模型在非淹 没网格中模拟污染物在陆地上的产生及随径流的输移过程；采用水体模型在 淹没网格中模拟污染物在水体中的对流扩散过程。
77.通过将陆地模型与水体模型直接动态双向耦合的方式，克服了采用单一 的陆地模型或水体模型以及陆地模型与水体模型单向耦合预测模拟流域水环 境污染的不足。通过根据预设水深阈值判断网格的状态，将水深较深的区域 视为作为淹没区，发生污染物的对流扩散过程，采用水体模型计算淹没区水 深和流速，并求解污染物的浓度，将水深较浅的区域作为非淹没区，发生污 染物的输移过程，采用陆地模型模拟污染物在陆地上的产生及输移过程，建 立的两个模型的动态双向耦合，计算中不指定边界位置，污染物和水流进入 淹没区的地点是动态变化的，洪水传播也可以淹没原有的汇流点，边界推移 形成新的耦合边界，实现了两种区域的交替转换，模拟更接近真实情境。对 比外部的单一耦合方式需要先后计算陆地模型和水体模型，本发明采用的双 向耦合方法中两种模型同时计算，污染物质可以由非淹没区汇入淹没区之中， 而河道湖泊中的污染物质也可以侵入汇水陆地区域，由此实现了岸上和水体 中的污染物的相互影响，更符合物理实际，更具科学性。耦合动边界作为控 制体单元的边界，在动边界两侧采用不同的计算模型(陆地模型和水体模型)， 通过水量、动量、污染物质量守恒的有限体积格式，既保证了各种物理量的 守恒，也可自动实现陆地污染物质通过耦合动边界进入水体，实现陆地模型 对水体模型的驱动，而不是通过人为地设置边界条件，将陆地污染物质作为 边界条件加载到水体模型中。
78.作为本发明的一个优选方案，根据预设水深阈值，对将要预测流域进行 网格划分，将水深低于预设水深阈值的网格作为非淹没网格，水深大于所述 预设水深阈值的网格作为淹没网格包括：
79.将将要预测流域划分为均匀的矩形网格；
80.对矩形网格的边界以及相邻的矩形网格之间的关系进行定义，并基于定 义对矩形网格进行编号；
81.将变量储存至编号的矩形网格的中心位置，给每个编号的矩形网格设置 初始条件，并将水深低于预设水深阈值的编号的矩形网格作为非淹没网格， 水深大于预设水深阈值的编号的矩形网格作为淹没网格。
82.采用均匀矩形网格，由于结构化网格具有实现容易、生成速度快、网格 质量好、数据结构简单的优点，可较好地与地形栅格数据相匹配，模型可直 接读取gis格式的数据信息，如dem、土地利用信息、土壤数据等，提高数 据前处理的效率，因此本发明的优选实施例采用结构化的均匀矩形网格离散 计算区域。首先将计算区域划分为均匀的矩形网格，设定预设水深阈值，根 据网格水深将计算域划分为淹没区和非淹没区，当网格的水深大于预先设定 的预设水深阈值，将该区域作为淹没区；若网格的水深小于该阈值，则将该 区域作为非淹没区。淹没区和非淹没区的范围随时间发生变化，两者的耦合 边界也处于动态变化中。不同区域采用不同的模型，在非淹没区采用陆地模 型模拟污染物的产生及输移过程；在淹没区应用水体模型模拟污染物在水中 的对流扩散过程。
83.作为本发明的一个优选方案，初始条件至少包括高程、水深、坡度、糙 率、外部污染负荷和污染物扩散系数，当然也可根据实际需要包括其它。
84.作为本发明的一个优选方案，陆地模型包括：
85.用于计算非淹没网格的汇流速度和水深的水文模型、用于根据非淹没网 格的汇流速度和水深，计算非淹没网格中的污染物的累积量和冲刷量的非点 源污染负荷模型和用于将非淹没网格中的污染物的累积量和冲刷量作为源 项，获得非点源污染物负荷分布数据的污染物输移模型；
86.水体模型包括：
87.用于将非淹没网格的汇流速度和水深作为水体水流运动的驱动条件，通 过求解局部区域的二维浅水方程，获得淹没网格内的流速和水深的水动力模 型和用于将非淹没网格的非点源污染物的浓度分布数据作为水体中污染物对 流扩散过程的驱动条件，并根据淹没网格内的流速和水深，求解污染物在水 体中的对流扩散过程的水质模型。
88.水文模型：水文模型包括产流和汇流过程。产流过程包括降雨、蒸发和 下渗，其中，即时计算中可忽略蒸发量，采用反距离插值法计算降雨量，基 于格林-安普顿(green-ampt)公式计算下渗量。产流过程计算的净雨量作为 质量守恒方程的源项。采用扩散波方程进行汇流计算。该模式数值稳定性好， 并能提高计算效率。其中，质量守恒方程为：
[0089][0090]
扩散波方程为：
[0091][0092][0093]
式中，h为水深，单位为m；q
x
和qy分别是x和y方向的网格单宽流量， 单位为m2·
s-1
；
qr为降雨量扣除入渗量的产流即净雨量，单位为m
·
s-1
；q
x
和qy分别为x和y方向的流量，单位为m3·
s-1
；a为过流断面的面积，单位为m2； r为水力半径，单位为m；s为水位坡度；n为地表糙率。
[0094]
非点源污染负荷模型：污染负荷模型包括污染物的累积和冲刷过程。干 旱无雨期间，污染物在土地上累积增长，采用指数函数描述污染物干旱时在 土地上的累积过程，认为污染物的累积遵从指数增长，直至最大值。土地上 累积的污染物在雨季时，将被雨水冲刷而进入径流中，进而进入管渠和河道， 直至土地上的污染物累积量被耗尽，采用指数冲刷函数计算污染物的冲刷量， 冲刷负荷正比于净流量的n次幂与污染物累积量的乘积。污染物累积方程公 式表示如下：
[0095][0096]
式中，bk为第k种污染物在单位面积的累积量，单位kg
·
ha-1
；c1为最大 累积量，单位kg
·
ha-1
，c2为累积速率常数，单位day-1
，参数c1,c2的选取可 以参考文献资料或通过试验确定，参数的取值随研究区域地点及污染物组分 不同而变化，c1的参考取值范围为1.28～12.44kg
·
ha-1
；c2的参考取值范围为 0.1～0.4day-1
；t是污染物累积时间，单位day，该累积方程表明污染物的增长 与时间成指数函数关系，当增长量达到最大时，增长停止。
[0097]
常规的非点源污染负荷量计算将在子流域(或子区域)的面积上进行积 分，如公式一所示，获得污染物在区域上的总量w随时间的变化，如公式二 所示，其中，公式一为：公式二为：
[0098]
式中，(bk)
total
为第k种污染物在流域上总的累积量，单位kg；ω为子区域 的面积，单位ha；w为冲刷负荷，单位kg/h；c3为冲刷系数，单位(mm
·
h-1
)-c4
·
h-1
； c4为冲刷指数，无量纲；q为单位面积的径流速率，单位mm
·
h-1
，的单 位为h-1
。
[0099]
采用(bk)
total
的陆地非点源污染负荷计算方法不能反映污染物在子区域的 空间及时间变化。为了反映非点源污染负荷随空间位置及时间的变化，本发 明计算单位面积的污染负荷，如公式三所示，其中，公式三为：
[0100][0101]
式中，s
nps,k
是第k种非点源污染物经降雨冲刷进入陆地，并随径流输移 的单位面积的负荷量，单位kg
·
m-2
·
s；其余符号的单位及意义同前。其中，参 数c3,c4的选取可以参考文献资料或通过试验确定，参数的取值随研究区域地 点及污染物组分不同而变化，c3的参考取值范围为2～9，c4的取值范围为 0.2～0.6。
[0102]
陆地上的污染物输移模型：采用水质方程模拟污染物在陆地上随径流的 输移过程，计算公式如公式四所示。其中，公式三计算得到的单位面积上的 污染物负荷作为该水质方程的源项。其中，公式四为：
[0103][0104]
式中，ck为第k种污染物的浓度，单位kg
·
m-3
，h为水深，单位m；u，v 分别为x和y方向的流速，单位m
·
s-1
，由水文模型计算得到；d
x,k
和d
y,k
分别 第k种污染物在x和y方向的扩散系数，单位m2·
s-1
；s
nps,k
为第k种污染物的 非点源负荷量，单位kg
·
m-2
·
s-1
；s
c,k
为第k种
污染物组分自身降解或者化学反 应生成的量，单位kg
·
m-3
·
s-1
，一般可以表示为s
c.k
＝-λhc，λ为水体的自净系 数。
[0105]
公式四通过增加陆地非点源源项，使得陆地污染物的空间分布计算成为 可能。在陆地模型与水体模型耦合时，在陆地上应用公式四计算陆地污染物 的输移过程后获得的污染物浓度信息，将作为水体中污染物对流扩散计算的 条件。这样的条件将会更加准确地刻画陆地污染物进入水体的实际空间位置， 为实现陆地模型与水体模型的动态双向耦合奠定了基础。
[0106]
水动力模型：水动力模型采用二维浅水方程(swe)计算，swe由连续 方程和动量方程组成。产流计算得到的水量变化作为连续方程的源项。应用 动量方程计算水流的运动过程。时间步长受库朗稳定性条件的限制。swe的 公式表示形式如下：
[0107][0108]
式中，u是守恒物理量，f、g分别是x和y方向的对流项，表述形式如 下：
[0109]
u＝[h,hu,hv]
t
，f＝[hu,hu2+gh2/2,huv]
t g＝[hv,huv,hv2+gh2/2]
t
[0110]sx
和sy分别是x和y方向的底坡源项和摩阻源项，表述形式如下：
[0111][0112]
式中，h为水深，单位m；u，v分别为x和y方向的流速，单位m
·
s-1
；z 是高程，单位m；c是表征地表糙率的谢才系数，由曼宁公式计算。
[0113]
水体中的水质模型：采用水质方程模拟污染物在水体中的对流扩散过程。 水动力模型计算得到的水深和流速作为已知量代入到水质方程中，进而求得 污染物的浓度。其中，陆地非点源污染物的冲刷量通过耦合动边界进入水体， 污染物自身的降解以及外部输入的污染物的通量作为方程的源项。忽略底泥 释放源项的水质方程作为公式五，其中，公式五如下：
[0114][0115]
式中，符号含义与公式四相同。与公式四相比，公式五忽略了底泥的污 染物释放量，即不考虑非点源源项s
nps,k
。
[0116]
由于公式四与公式五具有相同的形式，通过陆地模型与水体模型的内部 连接边界，进行陆地区域与水体区域的污染物的动态交换，保证在耦合动边 界上的水量及污染物成分的通量守恒。
[0117]
作为本发明的一个优选方案，在对于非淹没网格采用陆地模型模拟污染 物在陆地上的产生及随径流的输移过程，计算非淹没区的污染物浓度分布数 据的过程中，
[0118]
在水文模型中：
[0119]
通过反距离插值法计算非淹没网格流域的降雨量，通过格林-安普顿公式 计算流域的下渗量；
[0120]
根据流域的降雨量和下渗量计算非淹没网格的净雨量；
[0121]
采用扩散波方程计算非淹没网格的汇流速度和水深；
[0122]
在非点源污染负荷模型中，基于水文模型计算得到的汇流速度和水深， 采用指数方程计算随降雨径流进入水体的污染物的量；
[0123]
在污染物输移模型中：
[0124]
采用水质方程模拟污染物在陆地上随水流的输移过程，采用godunov中 心格式的有限体积算法对水质方程进行离散，得到非淹没区的污染物负荷分 布数据。
[0125]
作为本发明的一个优选方案，在对于非淹没网格采用陆地模型模拟污染 物在陆地上的产生及随径流的输移过程，计算非淹没区的污染物负荷分布数 据的过程中，
[0126]
在水文模型中：
[0127]
通过反距离插值法计算非淹没网格流域的降雨量，通过格林-安普顿公式 计算流域的下渗量；
[0128]
根据流域的降雨量和下渗量计算非淹没网格的净雨量；
[0129]
采用扩散波方程计算非淹没网格的汇流速度和水深；
[0130]
在非点源污染负荷模型中，基于水文模型计算得到的汇流速度和水深， 采用指数方程计算随降雨径流进入水体的污染物的量；
[0131]
在污染物输移模型中：
[0132]
采用水质方程模拟污染物在陆地上随水流的输移过程，采用godunov中 心格式的有限体积算法对水质方程进行离散，得到非淹没区的污染物负荷分 布数据。
[0133]
作为本发明的一个优选方案，在对于淹没网格采用水体模型模拟污染物 在水体中的对流扩散过程，并将非淹没区的污染物负荷分布数据作为计算条 件，计算淹没区的污染物的浓度分布数据的过程中，
[0134]
水动力模型与所述水质模型采用紧密耦合求解的方法，进行耦合求解， 采用muscl线性插值的方法将淹没网格中心点处的变量插值到淹没网格的 界面上；
[0135]
采用godunov格式的有限体积法进行空间离散，通过hllc近似黎曼求 解器求解淹没网格界面上的对流通量，采用二阶runge-kutta法保证时间积分 的二阶精度，采用muscl格式保证空间上的二阶精度，求解得到淹没区网 格界面的水深和流速；
[0136]
基于耦合动边界界面的质量、动量守恒，将非淹没区的污染物负荷分布 作为水体模型的计算条件，并将淹没网格界面的水深和流速作为已知量代入 到水质方程中，采用godunov格式进行空间离散，采用hllc近似黎曼解求 解通过淹没网格界面的污染物质的通量，采用二阶runge-kutta法保证时间积 分的二阶精度，采用muscl格式保证空间上的二阶精度，得到淹没区的污 染物的浓度分布数据。
[0137]
作为本发明的一个优选方案，模型计算的时间步长的确定方法包括：
[0138]
采用显式格式的网格离散方法，确定模型计算的时间步长；其中，时间 步长的限制条件为：
[0139][0140]
式中，c是常数，用于保证格式的稳定性，1/2≤c＜1；u
ni,j
,v
ni,j
分别为当 前单元在x和y方向的流速分量，单位m
·
s-1
；h
ni,j
为单元网格的水深，单位m； g为重力速度；d
x
和dy分别为x和y方向的扩散系数，单位m2·
s-1
；γ为扩散 项稳定性系数。
[0141]
网格离散方法采用显式格式进行，计算时间步长受cfl (courant-friefrichs-lewy)条件限制。在计算中，当采用不同属性的网格计 算时，为保证计算的同步性，在计算中采用统一的时间步长。
[0142]
作为本发明的一个优选方案，
[0143]
基于模型计算的时间步长，计算淹没网格和非淹没网格的水深数据，并 根据水深数据对所述耦合动边界进行更新，包括：
[0144]
计算淹没网格和非淹没网格间断处的特征波，并根据特征波的方向确定 淹没区和非淹没区之间的质量和动量传递方向；
[0145]
根据淹没区和非淹没区之间的质量和动量传递方向，更新非淹没区和淹 没区的水深和流速；
[0146]
根据非淹没区和淹没区的水深，将水深低于预设水深阈值的网格作为更 新的非淹没网格，水深大于预设水深阈值的网格作为更新的淹没网格；
[0147]
根据更新的非淹没网格所形成的新的非淹没区和更新的淹没网格所形成 的新的淹没区之间的边界更新连接所述陆地模型与所述水体模型的耦合动边 界。
[0148]
作为本发明的一个优选方案，计算淹没网格和非淹没网格间断处的特征 波，并根据特征波的方向确定淹没区和非淹没区之间的质量和动量传递方向 包括：
[0149]
当耦合动边界处的特征波的方向全部指向淹没网格，则根据非淹没网格 的水深和流速计算淹没网格与非淹没网格的间断处的通量，并判定非淹没区 向淹没区传递质量，不传递动量；
[0150]
当所述耦合动边界处的特征波在淹没网格和非淹没网格均有分布，则根 据淹没网格的水深及流速和非淹没网格的水深及流速计算淹没网格与非淹没 网格的间断处的通量，并判定非淹没区向所述淹没区传递质量，淹没区向非 淹没区传递动量；
[0151]
当耦合动边界处的特征波的方向全部指向非淹没网格，则根据淹没网格 的水深和流速计算淹没网格与非淹没网格的间断处的通量，并判定淹没区向 非淹没区传递质量和动量。
[0152]
其中，淹没区和非淹没区之间的质量传递通过水量平衡进行计算，动量 传递通过流速进行计算。
[0153]
本发明提供的一种流域水环境模拟装置，包括：
[0154]
网格划分模块，用于根据预设水深阈值，对将要预测流域进行网格划分， 将水深低于预设水深阈值的网格作为非淹没网格，水深大于预设水深阈值的 网格作为淹没网格；
[0155]
模型建立模块，用于以非淹没网格所形成的非淹没区与淹没网格所形成 的淹没区之间的边界为耦合动边界，建立陆地模型和水体模型；
[0156]
计算模块，用于对于非淹没网格采用陆地模型模拟污染物在陆地上的产 生及随径流的输移过程，计算非淹没区的污染物负荷分布数据；对于淹没网 格将非淹没区的污染物负荷分布数据作为计算条件，采用水体模型模拟污染 物在水体中的对流扩散过程，计算淹没区的污染物的浓度分布数据；
[0157]
网格更新模块，用于基于模型计算的时间步长，计算淹没网格和非淹没 网格的水深数据，并根据水深数据对耦合动边界进行更新。
[0158]
为了更好的对本发明提供的流域水环境模拟方法及装置进行说明，提供 的具体
实施例如下：
[0159]
实施例1
[0160]
如图2至图14共同所示，在模拟水环境之前，通过仪器测量或收集流域 水文资料，采集山区流域地形数据、土地利用类型数据、降雨过程数据以及 污染物计算相关数据。具体计算过程如下：
[0161]
1)网格离散及常规参数赋值。基于结构化网格，将计算区域划分为均匀 的矩形网格。图4定义了网格边界及相邻网格的关系，(i,j)为当前控制体编 号(或单元网格编号)，(i-1,j),(i+1,j)为当前单元网格在x方向的左、右(l， r)两侧相邻的单元网格，(i,j-1),(i,j+1)为当前单元y方向上面、下面(up， dn)单元。变量储存在单元中心位置。给每个网格设置初始条件，如高程、 水深、坡度、糙率、外部污染负荷、污染物扩散系数等。
[0162]
2)网格状态判断。设置预设水深阈值，当网格的水深小于阈值时，作为 非淹没区，采用陆地模型求解污染物在陆地上的输移过程；当网格的水深大 于阈值时，作为淹没区，应用水体模型模拟污染物在水体中的对流扩散过程。
[0163]
3)确定模型计算的时间步长δt。网格离散方法采用显式格式进行，计算 时间步长受cfl(courant-friefrichs-lewy)条件限制。在计算中，当采用不 同属性的网格计算时，为保证计算的同步性，在计算中采用统一的时间步长 δt。限制条件如下：
[0164][0165]
式中，c是常数，用于保持格式的稳定性，1/2≤c＜1；u
ni,j
,v
ni,j
分别当前单 元在x和y方向的流速分量，单位m
·
s-1
；h
ni,j
为单元格的水深，单位m；g为 重力速度；d
x
和dy分别为x和y方向的扩散系数，单位m2·
s-1
；γ为扩散项稳 定性系数。
[0166]
4)求解陆地模型
[0167]
①
水文模型计算。首先通过反距离插值计算得到流域的降雨量，按 green-ampt公式进行计算得到下渗量。汇流计算采用dwe，单元网格内水流 的流向为水位梯度最大的方向。采用有限差分法求解dwe：
[0168][0169][0170][0171]
式中，上角标n是时间步数，q
nx(i,j)
,q
ny(i,j)
分别x和y方向的陆地汇流流 量，单位m3·
s-1
；s
nx(i,j)
,s
ny(i,j)
分别x和y方向的水位坡度；h
ni,j
为网格(i,j)的有 效水深，单位m；为网格左右两侧界面的有效水深，单位m；指的是 网格上下界面的有效水深，单位m；为网格单元的水深，单位m， w为矩形网格边长(正方形网格时为同一值，m)；nr为曼宁系 数。单元格的示意图可参考图4。
[0172]
时间步长末，网格单元的水深按下式进行更新计算：
[0173][0174]
则单元内单个时间步长内的水深变化为：
[0175][0176]
式中，i为单元网格(控制体)编号；分别为单元网格(i,j)在当前 时刻和下一时刻的水深，单位m；a
i,j
为单元网格的面积；和分别 为网格(i,j)在当前时刻的入流、出流量以及源项，单位m3·
s-1
。
[0177]
②
污染物负荷计算。污染物负荷计算流程如图5所示。干旱无雨时，污 染物在地表累积；而当降雨来临时，这些富集在地表的污染物质在降雨径流 的冲刷下进入水体，从固结在地表的污染块变成溶解于水中的物质。基于水 文模型计算得到的汇流过程，采用指数方程描述污染物随降雨径流的进入水 体的污染物质的通量。
[0178]
③
模拟污染物在陆地上的输移过程。采用水质方程模拟污染物在陆地上 随水流的输移过程。采用godunov中心格式的有限体积算法对水质方程进行 离散，水质方程的离散式如下：
[0179][0180]
式中，i,j为控制体编号；m为当前控制体的边数，n为矩形网格的边数， 均匀的矩形网格情况下，n＝4；a
i,j
为(i,j)控制体的面积。其余符号及意义同 前。
[0181]
水质方程离散式对流通量计算方法可根据双曲方程的黎曼解(riemannsolver)，通过判断特征线的符号，计算通过界面的通污染物量。应用muscl 线性重构和二阶显式runge-kutta方法保证时空二阶精度。
[0182]
muscl插值格式如图6所示。变量存储在网格中心点处，而网格界面上 未存储变量，但计算的是网格界面上的对流通量，因此采用muscl线性插 值的方法将网格(i,j)和(i+1,j)中心点处的变量插值到界面(i+1/2,j)上。由图 6(c)所示，
[0183][0184][0185]
由此可见，u
i+1/2 l
≠u
i+1/2 r
[0186]
式中u表示储存在网格中心点所有变量，如水深，流速，污染物通量等。
[0187]
因此，在求解污染物输移过程时，采用muscl线性插值的方法将网格 中心的变量插值到界面(i+1/2,j)，得到：h
ni+1 2-,j
，z
ni+1/2-,j
，(ck)
ni+1/2 l,j
， h
ni+1/2+,j
，z
ni+1/2+,j
，(ck)
ni+1/2 r,j
和
[0188]
基于此，可以在网格界面上计算得到四条特征线，如图7所示。特征线 的计算方法
为界面(i,j+1/2)下侧的两条特征线，和为界面(i,j+1/2)上侧的两条 特征线，如图8所示。
[0212]
由此得到，界面(i,j+1/2)上的特征线为：
[0213][0214][0215][0216]
由此根据计算通过界面(i,j+1/2)的污染物通量为
[0217][0218]
由此得到通过界面的污染物通量，采用muscl重构的方法可保证在空 间上的二阶精度。
[0219]
利用水质方程的离散式求变量在(n+1)时刻数值时，对流项在时间方向 上采用两步计算，达到时间上的二阶精度。
[0220][0221][0222][0223]
式中，φ为对流项离散算子；w为待求变量。
[0224]
水质方程的离散式中的扩散项计算方法如下：
[0225][0226][0227]
求解陆地模型的污染物对流扩散方程中使用的水深和流速，由步骤
①
水 文模型计算得到。并将扩散项和步骤
②
计算得到的非点源污染物冲刷量，作 为源项加入到水质方程中，由此得到污染物在陆地上随水流的输移过程。
[0228]
5)求解水体模型。
[0229]
水体模型中的水动力方程和水质方程采用紧密耦合求解的方法，即在水 动力方
程和水质方程同时计算，将水动力方程的变量代入到水质方程中，进 行耦合求解，能更精确地反应水动力变化对水质等物理过程的影响，这也是 本发明的创新点之一。采用godunov中心格式的有限体积算法进行空间离散， 将实际连续的水流离散化后得到数值计算中不连续的间断分布，在每个时刻 求解间断的通量得到单元格的输入、输出通量，进而确定单元内部水流物理 量(水深、速度、浓度)随时间的变化。
[0230]
通过hllc近似黎曼求解器求解单元界面上的对流通量。首先采用 muscl线性插值的方法分别将将网格(i,j)和(i+1,j)中心点处的变量插值到 界面(i+1/2,j)上，得到和和y方向上的变量也采用该方法处理。
[0231]
守恒型式的二维水动力及水质方程耦合的矢量形式表示如下：
[0232][0233]
式中：
[0234][0235][0236]
为了获得保持静态平衡的良好数值格式，尤其是在静水中，采用了重构 技术，并对界面上的通量计算也进行了修改。因此在水体模型中，使用以下 三种方法：使用良好平衡方案(即静水压重构)求解底坡，使用半隐式方法 求解摩阻项，污染物扩散的求解与陆地模型相同。
[0237]
静水压力重构的目的是重构数值通量中使用的新变量，以保持静止状态 和水深的正性。静水压重构可以写成：
[0238]hi
±
1,j
＝max[0,h
i,j
+z
i,j-max(z
i,j
,zi±
1,j
)]
[0239]hi,j
±1＝max[0,h
i,j
+z
i,j-max(z
i,j
,z
i,j
±1)]
[0240]
基于静水重构技术，对通量计算进行了修正。在这种情况下，数值格式 中的通量梯度必须与源项的近似值平衡。为了分析这一点，只考虑源项中的 底坡，因为摩擦项不影响静水解。因此，首先计算了底坡和修正通量。二维 水动力及水质方程耦合形式可改写成如下：
[0241]
[0242]
式中，
[0243]
对该线性方程组进行离散，得到具有二阶时空精度的齐次解的计算公式:
[0244][0245]
其中：
[0246][0247][0248][0249][0250][0251][0252]
式中，h
i+1/2+,j
,h
i+1/2-,j
分别为界面(i+1/2,j)左右两侧的水深；h
i+1/2l,j
,h
i+1/2r,j
分别 为(i+1/2,j)界面左右两侧的有效水深，h
i,j+1/2 r
和h
i,j+1/2 l
分别为(i,j+1/2)界面上下 两侧的有效水深，计算方法同前；s为修正项，目的是为了保证静水平衡；sc 为中心值，目的是为了保证静水平衡的二阶精度。
[0253]
由于浅水流动的特殊性，f和g的四个变量采用不同的计算方法。 f(hu,hu2+gh2/2,huv,huc)的前2个变量(hu,hu2+gh2/2)和 g(hv,huv,hv2+gh2/2,hvc)的第1和第3个变量(hv，hv2+gh2/2)，可按hllc 近似黎曼求解器进行计算。
[0254]
x方向的通量计算如下：
[0255]
利用二维水动力模型计算得到的n时刻的流速和水深，分别在(i+1/2,j)界 面左右两侧得到4条特征线：和如图7所示。由此得到通过 界面(i+1/2,j)的通量如下：
[0256][0257]
同理，y方向的通量计算如下：
[0258]
利用二维水动力模型计算得到的n时刻的流速和水深，分别在(i,j+1/2)界 面上下两侧计算得到4条特征线：和如图8所示。由 此得到通过界面(i+1/2,j)的通量如下：
[0259][0260]
而f的第3个和第4个变量(huv,huck)和g的第2个和第4个变量 (huv,hvck)可直接由f和g的第一个变量(hu,hv)确定，计算方法如下：
[0261][0262][0263][0264]
[0265]
求变量在(n+1)时刻的数值时，对流项在时间方向上采用两步计算，达 到时间上的二阶精度。
[0266]
求解摩阻项和扩散项：
[0267]
采用半隐式的方法求解摩阻源项，有效的保证了算法的稳定性。
[0268][0269]
水质方程中的扩散项的处理方法与陆地模型一样。与陆地上的水质模型 不同，水体中的水质方程不计算非点源污染冲刷源项，只考虑由于水体自净 能力引起的污染物降解。
[0270]
6)根据获得的水深h更新下一时刻淹没区和非淹没区的计算范围，即当 h大于预设水深阈值，作为淹没区，采用水体模型计算，当h小于阈值，作为 非淹没区，采用陆地模型计算。由此重复步骤2)～5)，直至计算结束。
[0271]
如图9至图14，通过模拟二维坡地上的污染物运动过程来验证该模型 的性能，并验证坡道水文区域与水深较大的淹没区域的耦合作用。计算域长 为500米，宽为400米，其中x方向的坡度为s
x
＝0.02+0.0000149x，y方 向的坡度为s
x
＝0.02+0.0000116y。糙率其中 n
rx
＝0.1-0.0000168x，n
ry
＝0.1-0.0000168y。该算例的降雨随时间线性 变化，在t＝0s时，降雨强度为0mm/s，在t＝6000s时，降雨强度达到最大为 28.8mm/h，t＝12000s时，降雨强度又线性下降到0mm/s，总的模拟时长为 14400s。采用均匀的矩形网格离散该计算区域，如图10所示。初始条件为整 个区域内水深和流速为0，边界条件为x＝0和y＝0的两条边采用固壁边界，其 余两条边为自由出流边界。在该区域的右上角一圆形区域内设置初始污染条 件，圆心为(180，230)，半径为5m，如图11所示。污染物的质量为10g/m2， 给定扩散系数为0.5m2/s。
[0272]
首先将出口处的流量与其他模型的结果进行对比分析以验证该模型的合 理性。将结果与由jaber等开发的数值模型、以及中国水利水电科学研究院开 发的ifms开发的二维水动力模型进行了比较，结果如图12所示。该模型得 到的结果与水科院的ifms以及jaber等人的二维水动力模型结果相接近，表 明该发明的耦合模型可以较好的处理降雨产汇流以及水流在地形中的平面演 进问题。
[0273]
图13为本发明的流域水环境模型对变坡度变糙率算例进行产汇污模拟计 算的两个时刻污染物浓度分布图。该数值模拟体现了坡道上污染物质随降雨 径流的冲刷沿着坡道向下游运动，从图中还可以发现污染物质的扩散现象， 污染物质流体团的中心位置的浓度随着向流域下游的输移过程中，其浓度峰 值不断减小。坡道上游的水流运动的方向指向下游，但可以清楚的看到污染 物质向坡道横向扩散过程，扩散过程随着水流的紊动强度增大而增大，使污 染物质的分子运动越剧烈，所以在图13(b)中可以看到，污染物质浓度分布 更广。
[0274]
为了更好的说明本发明的流域水环境模型求解的守恒性质，统计了该算 例不同时刻留在流域内的污染物质量与排出流域污染物质量，如图14所示， 图中展示了本发明数
值计算具有良好的守恒性。随着降雨强度的加大，污染 物质流出流域主要在6000s～8000s之间，在6000s之前，降雨强度尚未达到峰 值，且污染物质流体团输移到出口位置也需时间。在8000s之后，留在流域内 部的污染物质较少，污染物质浓度较低。
[0275]
通过上述具体实施方式可看出，本发明提供的流域水环境模拟方法及装 置，通过将陆地模型与水体模型直接动态双向耦合的方式，克服了采用单一 的陆地模型或水体模型以及陆地模型与水体模型单向耦合预测流域水环境污 染的不足。通过根据预设水深阈值判断网格的状态，将水深较深的区域作为 淹没区，发生污染物的对流扩散过程，采用水体模型计算淹没区水深和流速， 并求解污染物的浓度，将水深较浅的区域作为非淹没区，发生污染物的输移 过程，采用陆地模型模拟污染物在陆地上的产生及输移过程，建立的两个模 型的动态双向耦合，计算中不指定边界位置，污染物和水流进入淹没区的地 点是动态变化的，洪水传播也可以淹没原有的汇流点，边界推移形成新的耦 合边界，实现了两种区域的交替转换，模拟更接近真实情境。对比外部的单 一耦合方式需要先后计算陆地模型和水体模型，本发明采用的双向耦合方法 中两种模型同时计算，污染物质可以由非淹没区汇入淹没区之中，而河道湖 泊中的污染物质也可以侵入汇水陆地区域，由此实现了岸上和水体中的污染 物的相互影响，更符合物理实际，更具科学性。耦合动边界作为控制体单元 的边界，在动边界两侧采用不同的计算模型(陆地模型和水体模型)，通过水 量、动量、污染物质量守恒的有限体积格式，既保证了各种物理量的守恒， 也可自动实现陆地污染物质通过耦合动边界进入水体，实现陆地模型对水体 模型的驱动，而不是通过人为地设置边界条件，将陆地污染物质作为边界条 件加载到水体模型中。
[0276]
如上参照附图以示例的方式描述了根据本发明提出的流域水环境模拟方 法及装置。但是，本领域技术人员应当理解，对于上述本发明所提出的流域 水环境模拟方法及装置，还可以在不脱离本发明内容的基础上做出各种改进。 因此，本发明的保护范围应当由所附的权利要求书的内容确定。
技术特征：1.一种流域水环境模拟方法，其特征在于，包括如下步骤：s1、根据预设水深阈值，对将要预测流域进行网格划分，将水深低于预设水深阈值的网格作为非淹没网格，水深大于所述预设水深阈值的网格作为淹没网格；s2、以所述非淹没网格所形成的非淹没区与所述淹没网格所形成的淹没区之间的边界为耦合动边界连接陆地模型和水体模型；s3、对于所述非淹没网格采用所述陆地模型模拟污染物在陆地上的产生及随径流的输移过程，计算所述非淹没区的污染物负荷分布数据；对于所述淹没网格将所述非淹没区的污染物负荷分布数据作为计算条件，采用所述水体模型模拟污染物在水体中的对流扩散过程，计算所述淹没区的污染物的浓度分布数据；s4、基于模型计算的时间步长，计算所述淹没网格和所述非淹没网格的水深数据，并根据所述水深数据对所述耦合动边界进行更新，返回步骤s2。2.根据权利要求1所述的流域水环境模拟方法，其特征在于，根据预设水深阈值，对将要预测流域进行网格划分，将水深低于预设水深阈值的网格作为非淹没网格，水深大于所述预设水深阈值的网格作为淹没网格包括：将所述将要预测流域划分为均匀的矩形网格；对所述矩形网格的边界以及相邻的所述矩形网格之间的关系进行定义，并基于所述定义对所述矩形网格进行编号；将变量储存至编号的矩形网格的中心位置，给每个编号的矩形网格设置初始条件，并将水深低于所述预设水深阈值的编号的矩形网格作为非淹没网格，水深大于所述预设水深阈值的编号的矩形网格作为淹没网格。3.根据权利要求1所述的流域水环境模拟方法，其特征在于，所述陆地模型包括：用于计算所述非淹没网格的汇流速度和水深的水文模型、用于根据所述非淹没网格的汇流速度和水深，计算所述非淹没网格中的污染物的累积量和冲刷量的非点源污染负荷模型和用于将所述非淹没网格中的污染物的累积量和冲刷量作为源项，获得非点源污染物负荷分布数据的污染物输移模型；所述水体模型包括：用于将所述非淹没网格的汇流速度和水深作为水体水流运动的驱动条件，通过求解局部区域的二维浅水方程，获得所述淹没网格内的流速和水深的水动力模型和用于将所述非淹没网格的非点源污染物的浓度分布数据作为水体中污染物对流扩散过程的驱动条件，并根据所述淹没网格内的流速和水深，求解污染物在水体中的对流扩散过程的水质模型。4.根据权利要求3所述的流域水环境模拟方法，其特征在于，在对于所述非淹没网格采用所述陆地模型模拟污染物在陆地上的产生及随径流的输移过程，计算所述非淹没区的污染物负荷分布数据的过程中，在所述水文模型中：通过反距离插值法计算所述非淹没网格流域的降雨量，通过格林-安普顿公式计算流域的下渗量；根据所述流域的降雨量和下渗量计算所述非淹没网格的净雨量；采用扩散波方程计算所述非淹没网格的汇流速度和水深；在所述非点源污染负荷模型中，基于所述水文模型计算得到的汇流速度和水深，采用
指数方程计算随降雨径流进入水体的污染物的量；在所述污染物输移模型中：采用水质方程模拟污染物在陆地上随水流的输移过程，采用godunov中心格式的有限体积算法对水质方程进行离散，得到所述非淹没区的污染物负荷分布数据。5.根据权利要求3所述的流域水环境模拟方法，其特征在于，在对于所述淹没网格采用所述水体模型模拟污染物在水体中的对流扩散过程，并将所述非淹没区的污染物负荷分布数据作为计算条件，计算所述淹没区的污染物的浓度分布数据的过程中，所述水动力模型与所述水质模型采用紧密耦合求解的方法，进行耦合求解，采用muscl线性插值的方法将所述淹没网格中心点处的变量插值到所述淹没网格的界面上；采用godunov格式的有限体积法进行空间离散，通过hllc近似黎曼求解器求解所述淹没网格界面上的对流通量，采用二阶runge-kutta法保证时间积分的二阶精度，采用muscl格式保证空间上的二阶精度，求解得到所述淹没区网格界面的水深和流速；基于所述耦合动边界界面的质量、动量守恒，将所述非淹没区的污染物负荷分布作为水体模型的计算条件，并将所述淹没网格的水深和流速作为已知量代入到水质方程中，采用godunov格式进行空间离散，采用hllc近似黎曼解求解通过所述淹没网格界面的污染物质的通量，采用二阶runge-kutta法保证时间积分的二阶精度，采用muscl格式保证空间上的二阶精度，得到所述淹没区的污染物的浓度分布数据。6.根据权利要求1所述的流域水环境模拟方法，其特征在于，所述模型计算的时间步长的确定方法包括：采用显式格式的网格离散方法，确定模型计算的时间步长；其中，所述时间步长的限制条件为：其中式中，c是常数，用于保证格式的稳定性，1/2≤c＜1；u
ni,j
,v
ni,j
分别为当前单元在x和y方向的流速分量(m
·
s-1
)；h
ni,j
为单元网格的水深(m)；g为重力速度；d
x
和d
y
分别为x和y方向的扩散系数(m2·
s-1
)；γ为扩散项稳定性系数。7.根据权利要求1所述的流域水环境模拟方法，其特征在于，基于模型计算的时间步长，计算所述淹没网格和所述非淹没网格的水深数据，并根据所述水深数据对所述耦合动边界进行更新，包括：计算所述淹没网格和所述非淹没网格间断处的特征波，并根据所述特征波的方向确定所述淹没区和所述非淹没区之间的质量和动量传递方向；根据所述淹没区和所述非淹没区之间的质量和动量传递方向，更新所述非淹没区和所述淹没区的水深和流速；根据所述非淹没区和所述淹没区的水深，将水深低于所述预设水深阈值的网格作为更新的非淹没网格，水深大于所述预设水深阈值的网格作为更新的淹没网格；根据所述更新的非淹没网格所形成的新的非淹没区和所述更新的淹没网格所形成的新的淹没区之间的边界更新连接所述陆地模型与所述水体模型的耦合动边界。8.根据权利要求7所述的流域水环境模拟方法，其特征在于，所述计算所述淹没网格和
所述非淹没网格间断处的特征波，并根据所述特征波的方向确定所述淹没区和所述非淹没区之间的质量和动量传递方向包括：当所述耦合动边界处的特征波的方向全部指向所述淹没网格，则根据所述非淹没网格的水深和流速计算所述淹没网格与所述非淹没网格的间断处的通量，并判定所述非淹没区向所述淹没区传递质量，不传递动量；当所述耦合动边界处的特征波在所述淹没网格和所述非淹没网格均有分布，则根据所述淹没网格的水深及流速和所述非淹没网格的水深及流速计算所述淹没网格与所述非淹没网格的间断处的通量，并判定所述非淹没区向所述淹没区传递质量，所述淹没区向所述非淹没区传递动量；当所述耦合动边界处的特征波的方向全部指向所述非淹没网格，则根据所述淹没网格的水深和流速计算所述淹没网格与所述非淹没网格的间断处的通量，并判定所述淹没区向所述非淹没区传递质量和动量。9.一种流域水环境模拟装置，其特征在于，包括：网格划分模块，用于根据预设水深阈值，对将要预测流域进行网格划分，将水深低于预设水深阈值的网格作为非淹没网格，水深大于所述预设水深阈值的网格作为淹没网格；模型建立模块，用于以所述非淹没网格所形成的非淹没区与所述淹没网格所形成的淹没区之间的边界为耦合动边界，建立陆地模型和水体模型；计算模块，用于对于所述非淹没网格采用所述陆地模型模拟污染物在陆地上的产生及随径流的输移过程，计算所述非淹没区的污染物负荷分布数据；对于所述淹没网格将所述非淹没区的污染物负荷分布数据作为计算条件，采用所述水体模型模拟污染物在水体中的对流扩散过程，计算所述淹没区的污染物的浓度分布数据；网格更新模块，用于基于模型计算的时间步长，计算所述淹没网格和所述非淹没网格的水深数据，并根据所述水深数据对所述耦合动边界进行更新。
技术总结本发明提供一种流域水环境模拟方法及装置，包括如下步骤：根据预设水深阈值，对将要预测流域进行网格划分，将水深低于预设水深阈值的网格作为非淹没网格，水深大于所述预设水深阈值的网格作为淹没网格；以非淹没网格与淹没网格之间的边界为耦合动边界连接陆地模型和水体模型；对于非淹没网格采用陆地模型模拟污染物在陆地上的产生及随径流的输移过程；对于淹没网格将非淹没区的污染物负荷分布数据作为计算条件，采用水体模型模拟污染物在水体中的对流扩散过程；根据水深数据对将要预测流域重新进行网格划分。利用本发明能够解决目前现有的水环境模拟方法割裂了流域水—陆水文、水质响应过程，不够科学，计算精度低等问题。计算精度低等问题。计算精度低等问题。
技术研发人员：申言霞 江春波 周琦 段艳华
受保护的技术使用者：清华大学
技术研发日：2022.01.27
技术公布日：2022/7/5
