基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法

摘要

本发明公开了一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法。构建湖泊三维水动力-水温-水质模型，将湖泊离散成若干网格单元，并对网格采用Arakawa？C模式布置变量，基于分裂算法将湖泊三维水动力-水温-水质模型中各算子按照其物理波动过程的波频快慢特性进行分类，对低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理，运用分裂算法分步离散求解模型，得到湖泊水域内不同位置和时间的三维流场、水温和水质指标浓度。本发明提出的基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法能较准确地反映湖泊水体中的动量迁移、热量传递和污染物输移等复杂的物理过程，具有计算稳定性好、计算精度和计算效率高的特点。

说明书

技术领域

本发明属于湖泊水环境数值模拟技术领域，更具体地，涉及一种基于分 裂算法的湖泊三维水动力-水温-水质模拟预测方法。本发明提出的湖泊三维 水动力-水温-水质模拟预测方法，可用于河湖连通工程、引水调度、湖泊水 体富营养化、湖泊水环境修复和管理等方面，为湖泊生态环境的治理和保 护提供有效的技术手段和科学的决策分析。

背景技术

湖泊作为陆地水资源的重要载体，是地球水圈的关键组成单元，是人类 赖以生存和可持续发展的自然依托，但目前我国湖泊却遭受着湖泊水资源 萎缩、水污染和富营养化等水环境问题，已给我国社会经济发展和人民的 生活环境带来了严重的影响。湖泊水环境保护和综合治理已经成为资源、 环境、生态领域亟需解决的关键科学问题。本发明所提出的湖泊三维水动 力-水温-水质模拟预测方法，对于促进我国水生态文明建设，打造天蓝、地 绿、水美的“美丽中国”和实现“中国梦”具有重要的科学意义。

原型观测和数值模拟是湖泊水环境研究的主要技术手段，但由于湖流运 动在时间和空间上呈现出高度复杂性，原型观测往往受条件限制，无法做 到同步监测，而数值模拟，由于费用低、物理意义明确，并能较真实地反 映湖泊水体的水动力、水温和水质状况和理化生过程，已经成为支持水环 境研究的重要技术手段。合适的湖泊水环境模型和稳定高效的模型求解算 法是湖泊水环境数值模拟技术的关键。目前，湖泊水环境数值模拟技术存 在以下不足：

1.传统湖泊水环境模型求解方法通常采用单一的数值格式对模型进行离 散求解，导致计算要求的稳定条件无法满足或计算工作量大、计算效率低 等问题。例如，采用纯显格式对模型进行离散，可能因不满足CFL稳定判 据而造成计算结果发散的现象；而采用纯隐格式对模型进行离散，虽然保 证了计算稳定性，但是必须在整个求解区域上解大规模非线性代数方程组， 导致计算效率低；

2.在实际问题中，模型的数学控制方程往往被简化成二维甚至是一维， 虽然合理的简化易于求解，但是自然水体是三维的，简化的一、二维模型 并不能真实地反映水体在纵向、横向、垂向上的变化过程；

3.湖泊水环境中物质输移、能量转换等物理过程联系紧密、相辅相成， 湖泊水动力过程是水温热量传递过程和污染物质输移过程的驱动条件。而 某些研究基于过于简单的假设，将水温、水质过程视为孤立的单一过程， 并未考虑它们与水动力等过程的耦合效应，导致模拟的结果与真实的湖体 情况存在较大差异。

随着科研人员对湖泊水体复杂过程更深刻的认识，国内外近年来发表的 与本发明相关的主要研究成果有：

《水力发电学报》2005年第1期《非静压假定的σ坐标下垂向二维浅水 模型的应用研究》对非静压假定的σ坐标下垂向二维浅水模型进行了应用 研究。该文采用的是垂向二维的浅水模型，而本发明构建的是三维水动力- 水温-水质模型。

《水动力学研究与进展A辑》2009年第24卷3期《基于非结构化网格 上的三维微幅表面流动非静压数值模型》建立了基于非结构化网格求解三 维微幅自由表面流动非静压数值模型，一种有限差分法和有限体积法相结 合的方法被用来在非结构化网格上离散控制方程。该文建立的非静压三维 模型是解决微幅自由表面的流动问题，而本发明解决的是湖泊水动力-水温- 水质模拟预测问题。

《环境科学学报》2012年第32卷12期《湖泊水质模型SALMO在太 湖梅梁湾的应用》对水质模型在太湖梅梁湾的运用做出了细致的研究。该 文的湖泊水质模型SALMO并没有考虑水动力等过程对水质过程的影响， 而本发明构建的水动力-水温-水质模型体现了水动力过程对水温及水质过 程的影响。

《昆明理工大学学报》杂志2013年38卷第一期《洱海湖泊及湖湾三维 水动力模型构建及特征分析》针对洱海湖泊及湖湾建立了三维水动力模型， 对其流场、水温进行了研究。虽然该文建立的模型适合描述湖泊三维水动 力过程，但是该文的模型求解方法并未涉及本发明的分裂算法思想。

《InternationalJournalforNumericalMethodsinFluids》杂志2008年56 卷第六期《Athree-dimensionalnon-hydrostaticverticalboundaryfittedmodel forfree-surfaceflows》建立了三维非静压有限体积模型，用来模拟在垂向边 界适应网格上的自由水面流动，而本发明是为了解决湖泊水动力-水温-水质 模拟预测问题。

中国专利公告号CN102156779A公告日是2011年8月17日名称为《地 下水流仿真与预测分析方法》，通过对采区地下水数据进行动态观测与采集， 利用数据引擎将地下水数据集成到图形工作站中，自动构建各水层动态水 位的有限元网格模型，同时确定参数分区及其相应的参数值，从而实现水 层水位动态模拟以及地下水运移仿真。该专利采用的是地下水模型，地下 水模型的物理机制不适用于本发明研究的湖泊水体。

综上所述，从目前国内外研究的发展趋势看，如何考虑湖泊水环境各物 理过程的耦合，如何解决湖泊水环境系统中各物理过程存在显著波频快慢 差异的问题，针对湖泊水环境复杂模型求解方法如何兼顾稳定性、计算精 度和计算效率的问题，这些都是亟待解决的技术难题。

发明内容

针对现有湖泊水环境模型亟需解决的上述技术问题和社会发展的需求， 本发明提出一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法， 基于分裂算法的思想，本发明将湖泊水环境系统中的复合物理波动过程视 作由若干不同波频的单一波动过程组成的集合，遵循地球物理问题一般将 低频波动过程和高频波动过程分类处理的准则，将湖泊三维水动力-水温- 水质模型中的各算子按照其物理波动过程的波频快慢特性，分为低频慢过 程算子和高频快过程算子两类，对其中的低频慢过程算子采用显式处理， 对高频快过程算子采用隐式处理，这种分裂算法解决了湖泊水环境系统中 因各物理波动过程存在显著波频差异，而难以选用合适的模型求解方法的 难题。本发明方法能对模型所有算子选用最适宜其物理特征的离散处理方 法，较为准确地反映湖泊水体中的动量迁移、热量传递和污染物输移等复 杂的物理过程，为湖泊三维水动力-水温-水质的数值模拟计算和预测提供了 有效的技术支持。

本发明提供了一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测 方法，包括以下步骤：

(1)构建湖泊三维水动力-水温-水质模型，模型由湖泊水动力控制方程 组、水温控制方程和水质控制方程组成，将湖泊计算域离散成若干网格单 元，并对网格采用ArakawaC模式布置变量，采用分裂算法首先将湖泊三 维水动力-水温-水质模型中的水动力控制方程组的各算子按照其物理波动 过程的波频快慢特性进行分类，对低频慢过程算子采用显式处理，对高频 快过程算子采用隐式处理，通过分裂算法分步离散并求解水动力控制方程 组，得到湖泊水域内不同位置和时间的三维流场u、v、w及其水深H，其 中u、v、w分别表示x、y、z方向的流速；

(2)通过步骤(1)求得湖泊三维流场u、v、w和水深H等水动力参数 的基础上，采用分裂算法分别将水温控制方程和水质控制方程中的算子分 为低频慢过程算子和高频快过程算子两类，对其中的低频慢过程算子采用 显式处理，对高频快过程算子采用隐式处理，通过分裂算法分步离散并求 解水温控制方程和水质控制方程，得到湖泊水域内不同位置和时间的水温T 和水质指标浓度C。

其中，步骤(1)包括以下子步骤：

(1.1)构建湖泊三维水动力-水温-水质模型，模型由反映湖泊水动力、 水温和水质复杂变化过程的一组数学物理控制方程组构成，包括：

水动力控制方程组：

连续方程${\partial }_{x}Hu+{\partial }_{y}Hv+{\partial }_{\sigma }\omega +{\partial }_t\eta =0---\left(1\right)$

水位方程${\partial }_t\eta +{\partial }_x\left(H{\int }_0^{1}ud\sigma \right)+{\partial }_y\left(H{\int }_0^{1}vd\sigma \right)=0---\left(2\right)$

x动量方程${\partial }_{t}u+u{\partial }_{x}u+v{\partial }_{y}u+\omega {\partial }_{\sigma }u=fv-g{\partial }_x\eta +{\partial }_x\left({\nu }_h{\partial }_{x}u\right)+\partial y\left({\nu }_h{\partial }_{y}u\right)+{\partial }_{\sigma }\left(\frac{{\nu }_v}{H^2}{\partial }_{\sigma }u\right)---\left(3\right)$

y动量方程${\partial }_{t}v+u{\partial }_{x}v+v{\partial }_{y}v+\omega {\partial }_{\sigma }v=-fu-g{\partial }_y\eta +{\partial }_x\left({\nu }_h{\partial }_{x}v\right)+{\partial }_y\left({\nu }_h{\partial }_{y}v\right)+{\partial }_{\sigma }\left(\frac{{\nu }_v}{H^2}{\partial }_{\sigma }v\right)---\left(4\right)$

水温控制方程：${\partial }_{t}T+u{\partial }_{x}T+v{\partial }_{y}T+\omega {\partial }_{\sigma }T={\partial }_x\left(A_{HT}{\partial }_{x}T\right)+{\partial }_y\left(A_{HT}{\partial }_{y}T\right)+{\partial }_{\sigma }\left(\frac{A_{VT}}{H^2}{\partial }_{\sigma }T\right)---\left(5\right)$

水质控制方程：${\partial }_{t}C+u{\partial }_{x}C+v{\partial }_{y}C+\omega {\partial }_{\sigma }C={\partial }_x\left(A_{HC}{\partial }_{x}C\right)+{\partial }_y\left(A_{VC}{\partial }_{y}C\right)+{\partial }_{\sigma }\left(\frac{v_b}{H^2}{\partial }_{\sigma }C\right)---\left(6\right)$

式中，$\sigma =\frac{z+h}{H}$是σ坐标变换；$\omega =\frac{w}{H}-\frac{\sigma }{H}\frac{\partial H}{\partial t}+\frac{u}{H}(\frac{\partial h}{\partial x}-\sigma \frac{\partial H}{\partial x})+\frac{v}{H}(\frac{\partial h}{\partial y}-\sigma \frac{\partial H}{\partial y})$是σ变换后的垂向流速；x，y为水平坐标，z为垂向坐标；h为自由水面到 水底的水深值，η为静水位偏移位移，H＝h+η为总水深；u，v分别为x，y 方向的流速；t为时间；T为水温；C为污染物浓度；Ω为地球 自转角速度，为所处纬度；V_h和V_v分别为水平涡黏系数和垂向涡粘系数； V_b为垂向扩散系数；g为重力加速度；A_HT和A_VT为水平温度扩散系数和垂 向温度扩散系数；A_HC和A_VC为水平水质扩散系数和垂向水质扩散系数；

(1.2)采用有限差分法离散模型控制方程组，对湖泊地形数据进行预处 理：利用小波包变换方法进行纹理特征提取，采用反距离加权平均插值法 进行空间插值获得湖底地形散点三维坐标；采用主成份分析方法对湖泊遥 感影像进行光谱特征提取，获取湖泊边界的平面二维坐标；

(1.3)根据得到的湖底地形散点三维坐标和湖泊边界的平面二维坐标， 将湖泊离散成若干网格单元，并对网格采用ArakawaC模式布置变量：对 一个正方体网格单元，规定i、j和k分别为x、y和z方向的网格编号索引， 变量ω设置在上下表面的中点上(i，j，k±1/2)，在每一个横截面中，变量 u设置在网格左右两边的中心(i±1/2，j，k)，变量v设置在前后两边的中 点上(i，j±1/2，k)，水深H、水质浓度C、水温T设置在网格中心(i，j， k)；

(1.4)采用分裂算法将水动力控制方程组中的水平动量方程各算子按照 其物理波动过程的波频快慢特性进行分类，缓行内重力波为低频慢过程， 表面重力长波为高频快过程，由此将水平动量方程(3)、(4)的算子分为 低频慢过程算子和高频快过程算子两类，其中，低频慢过程算子包括对流 项、科氏力项和水平涡粘项，高频快过程算子包括重力梯度项和垂向涡粘 项；

(1.5)对低频慢过程算子对流项、科氏力项和水平涡粘项做显式处理， 求得x方向中间流速u1和y方向中间流速v1，以x方向为例：

$\begin{array}{c}u{1}_{i+1/2,j,k}=-\frac{\Delta t}{2\Delta x}{u}_{i+1/2,j,k}^n\left(u_{i+3/2,j,k}^n-u_{i-1/2,j,k}^n\right)-\frac{v_{i+1/2,j,k}^{n*}\Delta t}{2\Delta y}\left(u_{i+1/2,j+1,k}^n-u_{i+1/2,j-1,k}^n\right)\\ -\frac{{\omega }_{i+1/2,j,k}^{n*}\Delta t}{{\mathrm{\Delta \sigma }}_k+0.5\left({\mathrm{\Delta \sigma }}_{k+1}+{\mathrm{\Delta \sigma }}_{k-1}\right)}\left(u_{i+1/2,j,k+1}^n-u_{i+1/2,j,k-1}^n\right)+\frac{v_h\Delta t}{{\mathrm{\Delta x}}^2}\left(u_{i+3/2,j,k}^n-u_{i-1/2,j,k}^n\right)\\ +\frac{v_h\Delta t}{{\mathrm{\Delta y}}^2}\left(u_{i,j+1,k}^n-2u_{i,j,k}^n+u_{i,j-1,k}^n\right)+\left(1-\frac{2v_h\Delta t}{{\mathrm{\Delta x}}^2}\right)u_{i+1/2,j,k}^n-{\mathrm{fv}}_{i+1/2,j,k}^{n*}\end{array}---\left(7\right)$

式中：

$v_{i+1/2,j,k}^{n*}=\frac{1}{4}(v_{i,j-1/2,k}^n+v_{i,j+1/2,k}^n+v_{i+1,j-1/2,k}^n+v_{i+1,j+1/2,k}^n)---\left(8\right)$

${\omega }_{i+1/2,j,k}^{n*}=\frac{1}{4}({\omega }_{i,j,k+1/2}^n+{\omega }_{i,j,k-1/2}^n+{\omega }_{i+1,j,k+1/2}^n+{\omega }_{i+1,j,k-1/2}^n)---\left(9\right)$

(1.6)由显式离散水位方程(2)，求得下一个时间节点的静水位偏移位 移ηⁿ⁺¹，根据已经求到的中间流速u1、v1和静水位偏移位移ηⁿ⁺¹，对高频 快过程算子重力梯度项做隐式处理，求得x方向中间流速u2和y方向中间 流速v2，以x方向为例：

$\frac{u{2}_{i+1/2,j,k}-u{1}_{i+1/2,j,k}}{\Delta t}=-g\frac{{\eta }_{i+1,j}^{n+1}-{\eta }_{i,j}^{n+1}}{\Delta x}---\left(10\right)$

(1.7)根据已经求到的中间流速u2和v2，对高频快过程算子垂向涡粘 项做隐式处理，整理后分别得到关于下一个时间节点的x方向流速uⁿ⁺¹和y 方向流速vⁿ⁺¹的线性方程式，以x方向为例：

$UT·u_{k+1}^{n+1}+UB·u_{k-1}^{n+1}+UC·u_k^{n+1}=UF---\left(11\right)$

式中UT、UB、UC、UF是包含已知变量u2的参数，在所有网格上由 式(11)得到一个涉及三维流场时空关系的线性矩阵方程组，采用托马斯

法求解矩阵方程组，得到uⁿ⁺¹，并用同样的方法求出vⁿ⁺¹；

(1.8)根据已经得到的x方向流速uⁿ⁺¹、y方向流速vⁿ⁺¹和连续方程(1)， 求得垂向流速ωⁿ⁺¹，根据σ反坐标变换求得笛卡尔坐标系下垂向流速wⁿ⁺¹；

(1.9)重复步骤(1.5)-(1.8)，循环迭代计算，得到湖泊水域内不同 位置和时间的三维流场u、v、w及其静水位偏移位移η，并根据静水位偏移 位移η进一步求得水深H。

其中，步骤(2)包括以下子步骤：

(2.1)根据步骤(1.4)类似的原理，将水温控制方程(5)和水质控制 方程(6)的算子分为低频慢过程算子和高频快过程算子两类，其中，低频 慢过程算子包括对流项、水平扩散项，高频快过程算子包括垂向扩散项；

(2.2)对低频慢过程算子对流项和水平扩散项做显式处理，求得中间温 度T1和中间浓度C1；

(2.3)根据求得的T1、C1，对高频快过程算子垂向扩散项做隐式处理， 求得下一个时间节点的温度Tⁿ⁺¹、浓度Cⁿ⁺¹，整理得到一个线性方程式，以 温度T为例：

${\mathrm{MT}}_{i,j,k}^{n+1}+{\mathrm{ST}}_{i,j,k+1}^{n+1}+{\mathrm{PT}}_{i,j,k-1}^{n+1}=R---\left(12\right)$

M、S、P、R是包含已知变量T1的参数，在各个网格写出式(12)，分 别得到一个涉及水温时空关系和水质时空关系的线性矩阵方程组，采用托 马斯法求得Tⁿ⁺¹，并用同样的方法求得Cⁿ⁺¹；

(2.4)重复步骤(2.2)-(2.3)，循环迭代计算，得到湖泊水域内不同 位置和时间的水温T和水质指标浓度C。

总体而言，本发明与现有技术相比，具有以下特点：

(1)基于分裂算法的思想，本发明将湖泊水环境系统中的复合物理波 动过程视作由若干不同波频的单一波动过程组成的集合，遵循地球物理问 题一般将低频波动过程和高频波动过程分类处理的准则，将湖泊三维水动 力-水温-水质模型的各算子按照其物理波动过程的波频快慢特性，分为低频 慢过程算子和高频快过程算子两类，对其中的低频慢过程算子采用显式处 理，对高频快过程算子采用隐式处理；

(2)能对模型所有算子选用最适宜其物理特征的离散处理方式，解决 了湖泊水环境系统中因各物理波动过程存在显著波频差异，而难以选用合 适的模型求解方法的技术难题；

(3)能较准确地反映湖泊水体中的动量迁移、热量传递和污染物输移 等复杂的物理过程，具有计算稳定性好、计算精度和计算效率高的特点。

附图说明

图1是本发明基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法 的总流程图；

图2是湖泊网格化及前处理过程图；

图3是本发明实施例中沙湖边界线示意图；

图4是本发明实施例中沙湖网格划分的示意图；

图5是本发明实施例中ArakawaC变量布置示意图；

图6是本发明实施例中湖泊某时刻的三维流场矢量图；

图7是本发明实施例中某网格水温时间序列图(垂向分为4层)；

图8是本发明实施例中湖泊某时刻COD浓度分布渲染图(垂向分为4 层)；

图9是本发明实例中静水位偏移位移η的实测值与模拟值；

图10是本发明实例中水温T的实测值与模拟值。

具体实施方式

本发明以武汉市沙湖为例。沙湖是武汉市重要的城市湖泊，分为内沙湖 (0.134平方公里)和外沙湖(3.197平方公里)。根据遥感影像对沙湖进行 特征提取，构建沙湖三维水动力-水温-水质模型。

图1所示为本发明一种基于分裂算法的湖泊三维水动力-水温-水质模拟 预测方法的总流程图，其总体思路是对湖泊水动力-水温-水质模型中的各算 子根据其波频特性进行分类离散处理。包括以下步骤：

(1)构建湖泊三维水动力-水温-水质模型，模型由湖泊水动力控制方程 组、水温控制方程和水质控制方程组成，将湖泊计算域离散成若干网格单 元，并对网格采用ArakawaC模式布置变量，采用分裂算法首先将湖泊三 维水动力-水温-水质模型中的水动力控制方程组的各算子按照其物理波动 过程的波频快慢特性进行分类，对低频慢过程算子采用显式处理，对高频 快过程算子采用隐式处理，通过分裂算法分步离散并求解水动力控制方程 组，得到湖泊水域内不同位置和时间的三维流场u、v、w及其水深H，其 中u、v、w分别表示x、y、z方向的流速；

具体而言，本步骤包括以下子步骤：

(1.1)构建湖泊三维水动力-水温-水质模型，模型由反映湖泊水动力、 水温和水质复杂变化过程的一组数学物理控制方程组构成，包括：

水动力控制方程组：

连续方程${\partial }_{x}Hu+{\partial }_{y}Hv+{\partial }_{\sigma }\omega +{\partial }_t\eta =0---\left(1\right)$

水位方程${\partial }_t\eta +{\partial }_x\left(H{\int }_0^{1}ud\sigma \right)+{\partial }_y\left(H{\int }_0^{1}vd\sigma \right)=0---\left(2\right)$

x动量方程${\partial }_{t}u+u{\partial }_{x}u+v{\partial }_{y}u+\omega {\partial }_{\sigma }u=fv-g{\partial }_x\eta +{\partial }_x\left({\nu }_h{\partial }_{x}u\right)+\partial y\left({\nu }_h{\partial }_{y}u\right)+{\partial }_{\sigma }\left(\frac{{\nu }_v}{H^2}{\partial }_{\sigma }u\right)---\left(3\right)$

y动量方程${\partial }_{t}v+u{\partial }_{x}v+v{\partial }_{y}v+\omega {\partial }_{\sigma }v=-fu-g{\partial }_y\eta +{\partial }_x\left({\nu }_h{\partial }_{x}v\right)+{\partial }_y\left({\nu }_h{\partial }_{y}v\right)+{\partial }_{\sigma }\left(\frac{{\nu }_v}{H^2}{\partial }_{\sigma }v\right)---\left(4\right)$

水温控制方程：${\partial }_{t}T+u{\partial }_{x}T+v{\partial }_{y}T+\omega {\partial }_{\sigma }T={\partial }_x\left(A_{HT}{\partial }_{x}T\right)+{\partial }_y\left(A_{HT}{\partial }_{y}T\right)+{\partial }_{\sigma }\left(\frac{A_{VT}}{H^2}{\partial }_{\sigma }T\right)---\left(5\right)$

水质控制方程：${\partial }_{t}C+u{\partial }_{x}C+v{\partial }_{y}C+\omega {\partial }_{\sigma }C={\partial }_x\left(A_{HC}{\partial }_{x}C\right)+{\partial }_y\left(A_{VC}{\partial }_{y}C\right)+{\partial }_{\sigma }\left(\frac{v_b}{H^2}{\partial }_{\sigma }C\right)---\left(6\right)$

式中，$\sigma =\frac{z+h}{H}$是σ坐标变换；$\omega =\frac{w}{H}-\frac{\sigma }{H}\frac{\partial H}{\partial t}+\frac{u}{H}(\frac{\partial h}{\partial x}-\sigma \frac{\partial H}{\partial x})+\frac{v}{H}(\frac{\partial h}{\partial y}-\sigma \frac{\partial H}{\partial y})$是σ变换后的垂向流速；x，y为水平坐标，z为垂向坐标；h为自由水面到 水底的水深值，η为静水位偏移位移，H＝h+η为总水深；u，v分别为x，y 方向的流速；t为时间；T为水温；C为污染物浓度；Ω为地球 自转角速度，为所处纬度；V_h和V_v分别为水平涡黏系数和垂向涡粘系数； V_b为垂向扩散系数；g为重力加速度；A_HT和A_VT为水平温度扩散系数和垂 向温度扩散系数；A_HC和A_VC为水平水质扩散系数和垂向水质扩散系数；

在本实施方案中，沙湖处于北纬29°58′，东经114°33′；水平涡黏 系数V_h和垂向涡粘系数V_v分别为1.5×10^-6和1.0×10^-4；水平水温分子扩散 系数A_HT和垂向温度扩散系数A_VT皆为6×10^-6；水平水质分子扩散系数A_HC和垂向水质扩散系数A_VC皆为2×10^-5；距水面10m处风速为2-5m/s；湖底 摩擦系数取0.025；重力加速度g＝10m/s²；水体比热容c_p＝4.2×10³J/(kg℃)；

(1.2)如图2所示，采用有限差分法离散模型控制方程组，对沙湖地形 数据进行预处理：利用小波包变换方法进行纹理特征提取，采用反距离加 权平均插值法进行空间插值获得沙湖湖底地形散点三维坐标，然后将沙湖 垂向分层为四层；

(1.3)通过谷歌地图获取沙湖的遥感影像，采用主成份分析方法对遥感 影像进行光谱特征提取，获取沙湖边界的平面二维坐标，生成湖泊边界线L (图3)；根据湖泊边界线L的4个极点坐标，生成湖泊边界线L的一个外 接矩形，将外接矩形划分为若干个方形网格。以网格中心点为标识，将网 格中心与湖泊边界线L进行空间拓扑分析，中心点在湖泊边界线L内的网 格保留，反之则删除，最终的网格数为305个(图4)；

(1.4)在划分好的网格上采用ArakawaC模式布置模型变量(图5)： 对一个正方体网格单元，规定i、j和k分别为x、y和z方向的网格编号索 引，变量ω设置在上下表面的中点上(i，j，k±1/2)，在每一个横截面中， 变量u设置在网格左右两边的中心(i±1/2，j，k)，变量v设置在前后两边 的中点上(i，j±1/2，k)，水深H、水质浓度C、水温T设置在网格中心(i， j，k)；

(1.5)采用分裂算法将水动力控制方程组中的水平动量方程各算子按照 其物理波动过程的波频快慢特性进行分类，缓行内重力波为低频慢过程， 表面重力长波为高频快过程，由此将水平动量方程(3)、(4)的算子分为 低频慢过程算子和高频快过程算子两类，其中，低频慢过程算子包括对流 项、科氏力项和水平涡粘项，高频快过程算子包括重力梯度项和垂向涡粘 项；

(1.6)对低频慢过程算子对流项、科氏力项和水平涡粘项做显式处理， 求得x方向中间流速u1和y方向中间流速v1，以x方向为例：

$\begin{array}{c}u{1}_{i+1/2,j,k}=-\frac{\Delta t}{2\Delta x}{u}_{i+1/2,j,k}^n\left(u_{i+3/2,j,k}^n-u_{i-1/2,j,k}^n\right)-\frac{v_{i+1/2,j,k}^{n*}\Delta t}{2\Delta y}\left(u_{i+1/2,j+1,k}^n-u_{i+1/2,j-1,k}^n\right)\\ -\frac{{\omega }_{i+1/2,j,k}^{n*}\Delta t}{{\mathrm{\Delta \sigma }}_k+0.5\left({\mathrm{\Delta \sigma }}_{k+1}+{\mathrm{\Delta \sigma }}_{k-1}\right)}\left(u_{i+1/2,j,k+1}^n-u_{i+1/2,j,k-1}^n\right)+\frac{v_h\Delta t}{{\mathrm{\Delta x}}^2}\left(u_{i+3/2,j,k}^n-u_{i-1/2,j,k}^n\right)\\ +\frac{v_h\Delta t}{{\mathrm{\Delta y}}^2}\left(u_{i,j+1,k}^n-2u_{i,j,k}^n+u_{i,j-1,k}^n\right)+\left(1-\frac{2v_h\Delta t}{{\mathrm{\Delta x}}^2}\right)u_{i+1/2,j,k}^n-{\mathrm{fv}}_{i+1/2,j,k}^{n*}\end{array}---\left(7\right)$

式中：

$v_{i+1/2,j,k}^{n*}=\frac{1}{4}(v_{i,j-1/2,k}^n+v_{i,j+1/2,k}^n+v_{i+1,j-1/2,k}^n+v_{i+1,j+1/2,k}^n)---\left(8\right)$

${\omega }_{i+1/2,j,k}^{n*}=\frac{1}{4}({\omega }_{i,j,k+1/2}^n+{\omega }_{i,j,k-1/2}^n+{\omega }_{i+1,j,k+1/2}^n+{\omega }_{i+1,j,k-1/2}^n)---\left(9\right)$

(1.7)由显式离散水位方程(2)，求得下一个时间节点的静水位偏移位 移ηⁿ⁺¹，根据已经求到的中间流速u1、v1和静水位偏移位移ηⁿ⁺¹，对高频 快过程算子重力梯度项做隐式处理，求得x方向中间流速u2和y方向中间 流速v2，以x方向为例：

$\frac{u{2}_{i+1/2,j,k}-u{1}_{i+1/2,j,k}}{\Delta t}=-g\frac{{\eta }_{i+1,j}^{n+1}-{\eta }_{i,j}^{n+1}}{\Delta x}---\left(10\right)$

(1.8)根据已经求到的中间流速u2和v2，对高频快过程算子垂向涡粘 项做隐式处理，整理后分别得到关于下一个时间节点的x方向流速uⁿ⁺¹和y 方向流速vⁿ⁺¹的线性方程式，以x方向为例：

$UT·u_{k+1}^{n+1}+UB·u_{k-1}^{n+1}+UC·u_k^{n+1}=UF---\left(11\right)$

式中UT、UB、UC、UF是包含已知变量u2的参数，在所有网格上由式 (11)得到一个涉及三维流场时空关系的线性矩阵方程组，采用托马斯法 求解矩阵方程组，得到uⁿ⁺¹，并用同样的方法求出vⁿ⁺¹；

(1.9)根据已经得到的x方向流速uⁿ⁺¹、y方向流速vⁿ⁺¹和连续方程(1)， 求得垂向流速ωⁿ⁺¹，根据σ反坐标变换求得笛卡尔坐标系下垂向流速wⁿ⁺¹；

(1.10)重复步骤(1.6)-(1.9)，循环迭代计算，得到湖泊水域内不同 位置和时间的三维流场u、v、w及其静水位偏移位移η，并根据静水位偏移 位移η进一步求得水深H。

(2)通过步骤(1)求得湖泊三维流场u、v、w和水深H等水动力参数 的基础上，采用分裂算法分别将水温控制方程和水质控制方程中的算子分 为低频慢过程算子和高频快过程算子两类，对其中的低频慢过程算子采用 显式处理，对高频快过程算子采用隐式处理，通过分裂算法分步离散并求 解水温控制方程和水质控制方程，得到湖泊水域内不同位置和时间的水温T 和水质指标浓度C。

具体而言，本步骤包括以下子步骤：

(2.1)根据步骤(1.5)类似的原理，将水温控制方程(5)和水质控制 方程(6)的算子分为低频慢过程算子和高频快过程算子两类，其中，低频 慢过程算子包括对流项、水平扩散项，高频快过程算子包括垂向扩散项；

(2.2)对低频慢过程算子对流项和水平扩散项做显式处理，求得到中间 温度T1和中间浓度C1；

(2.3)根据求得的T1、C1，对高频快过程算子垂向扩散项做隐式处理， 求得下一个时间节点的温度Tⁿ⁺¹、浓度Cⁿ⁺¹，整理得到一个线性方程式，以 温度T为例：

${\mathrm{MT}}_{i,j,k}^{n+1}+{\mathrm{ST}}_{i,j,k+1}^{n+1}+{\mathrm{PT}}_{i,j,k-1}^{n+1}=R---\left(12\right)$

M、S、P、R是包含已知变量T1的参数，在各个网格写出式(12)，分 别得到一个涉及水温时空关系和水质时空关系的线性矩阵方程组，采用托 马斯法求得Tⁿ⁺¹，并用同样的方法求得Cⁿ⁺¹；

(2.4)重复步骤(2.2)～(2.3)，循环迭代计算，得到沙湖内不同位置 和时间的水温T和水质指标浓度C；

(2.5)根据求得的三维流场u、v、w、水深H、水温T和水质指标浓度 C进行可视化后处理，具体为：

(a)根据某时刻的湖区三维流场的计算结果生成流场矢量图(图6)；

(b)对某个具体的网格，以时间为横坐标，以湖泊水环境某因子(流 场、水温、水质指标)的值为纵坐标，生成该因子的时间序列图(图7)；

(c)根据某时刻的湖区垂向每层的水温(或水质)计算结果生成三维 水温(或水质)分布渲染图(图8)。

模型验证

计算值与实测值对比

图9、图10分别是静水位偏移位移η、水温T的模拟值与实际观测数据 的对比结果；表1是水质指标(COD)浓度C的模拟值与实际观测数据的 对比结果：

表1沙湖水质指标COD模拟计算结果

计算稳定性与计算效率

本案例所有算法程序均运行于同一台PC(CPU:Celeron(R)E3300 2.50GHz；Memory:2GB,DDR2)。IDE(编译器)为visualstudio2010(c#)。 采用本发明方法与传统的纯显格式处理方法对模型进行求解，计算结果表 明纯显格式不能满足稳定条件，出现了计算结果发散的现象，而本发明方 法在整个计算过程中保持着良好的稳定性。在网格数(305)和时间步数(100) 的情况下，分别采用本发明方法与传统的纯隐格式处理方法对模型进行求 解，本发明方法和纯隐格式的CPU耗时分别为1.3和3.6个机时单位，本发 明方法比纯隐格式的计算效率提高了64％。

本发明方法的模拟结果与实测值吻合良好，应用实例证明本发明方法具 有良好的稳定性，且计算精度和计算效率较高。综上所述，本发明提出的 一种基于分裂算法的湖泊三维水动力-水温-水质模拟预测方法能较为准确 地反映湖泊水体中的动量迁移、热量传递和污染物输移等复杂的物理过程， 具有计算稳定性好、计算精度和计算效率较高的特点。