计算域离散成若干网格单元，并对网格采用Arakawa C模式布置变量，采用分裂算法首先将湖泊三维水动力-水温-水质模型中的水动力控制方 程组的各算子按照其物理波动过程的波频快慢特性进行分类，对低频慢过程算子采用显式 处理，对高频快过程算子采用隐式处理，通过分裂算法分步离散并求解水动力控制方程组， 得到湖泊水域内不同位置和时间的三维流场u、v、w及其水深H，其中u、v、w分别表示x、y、 z方向的流速；
[0069] 具体而言，本步骤包括以下子步骤：
[0070] (I. 1)构建湖泊三维水动力-水温-水质模型，模型由反映湖泊水动力、水温和水 质复杂变化过程的一组数学物理控制方程组构成，包括：
[0071] 水动力控制方程组：
[0072]
[0073]
[0074]
[0075]
[0076]
[0077]
[0078] 是σ变换后的垂向流速；X，y为水平坐标，Z为垂向坐标；h为自由水面到水底的水深值， η为静水位偏移位移，H = h+ η为总水深；U，V分别为X，y方向的流速；t为时间；T为水 温；C为污染物浓度；/=ZiJsinw Ω为地球自转角速度，於为所处炜度；Vjp V 别为水 平涡黏系数和垂向涡粘系数；vb为垂向扩散系数；g为重力加速度；A H#P A ντ为水平温度扩 散系数和垂向温度扩散系数;AHjP AveS水平水质扩散系数和垂向水质扩散系数；
[0079] 在本实施方案中，沙湖处于北炜29° 5V，东经114° 33';水平涡黏系数VjP 垂向涡粘系数Vv分别为1.5X10 6和1.0X10 4;水平水温分子扩散系数Aht和垂向温度扩 散系数Avt皆为6X10 6;水平水质分子扩散系数A Ηε和垂向水质扩散系数Ave皆为2X10 5; 距水面IOm处风速为2-5m/s ;湖底摩擦系数取0. 025 ;重力加速度g = 10m/s2;水体比热容 Cp= 4. 2X103J/(kg〇C )；
[0080] (1. 2)如图2所示，采用有限差分法离散模型控制方程组，对沙湖地形数据进行预 处理：利用小波包变换方法进行纹理特征提取，采用反距离加权平均插值法进行空间插值 获得沙湖湖底地形散点三维坐标，然后将沙湖垂向分层为四层；
[0081] (1. 3)通过谷歌地图获取沙湖的遥感影像，采用主成份分析方法对遥感影像进行 光谱特征提取，获取沙湖边界的平面二维坐标，生成湖泊边界线L (图3);根据湖泊边界线L 的4个极点坐标，生成湖泊边界线L的一个外接矩形，将外接矩形划分为若干个方形网格。 以网格中心点为标识，将网格中心与湖泊边界线L进行空间拓扑分析，中心点在湖泊边界 线L内的网格保留，反之则删除，最终的网格数为305个（图4);
[0082] (1. 4)在划分好的网格上米用Arakawa C模式布置模型变量（图5):对一个正方 体网格单元，规定i、j和k分别为X、y和z方向的网格编号索引，变量ω设置在上下表面 的中点上（i，j，k±l/2)，在每一个横截面中，变量u设置在网格左右两边的中心（i±l/2， j，k)，变量V设置在前后两边的中点上（i，j ± 1/2, k)，水深Η、水质浓度C、水温T设置在网 格中心（i，j，k);
[0083] (I. 5)采用分裂算法将水动力控制方程组中的水平动量方程各算子按照其物理波 动过程的波频快慢特性进行分类，缓行内重力波为低频慢过程，表面重力长波为高频快过 程，由此将水平动量方程（3)、⑷的算子分为低频慢过程算子和高频快过程算子两类，其 中，低频慢过程算子包括对流项、科氏力项和水平涡粘项，高频快过程算子包括重力梯度项 和垂向祸粘项；
[0084] (1. 6)对低频慢过程算子对流项、科氏力项和水平涡粘项做显式处理，求得X方向 中间流速ul和y方向中间流速vl，以X方向为例：
[0085]
[0086] 式中：
[0087] UiN 丄 丄 OUlO乙 Λ J ^ 〇/ ?
[0088]
[0089] (I. 7)由显式离散水位方程（2)，求得下一个时间节点的静水位偏移位移iT+1，根 据已经求到的中间流速ul、Vl和静水位偏移位移ιΓ +1，对高频快过程算子重力梯度项做隐 式处理，求得X方向中间流速u2和y方向中间流速ν2,以X方向为例：
[0090]
[0091] (1. 8)根据已经求到的中间流速U2和v2,对高频快过程算子垂向涡粘项做隐式处 理，整理后分别得到关于下一个时间节点的X方向流速u n+1和y方向流速V n+1的线性方程 式，以X方向为例：
[0092]
[0093] 式中UT、UB、UC、UF是包含已知变量u2的参数，在所有网格上由式（11)得到一个 涉及三维流场时空关系的线性矩阵方程组，采用托马斯法求解矩阵方程组，得到ιΓ 1，并用 同样的方法求出，+1;
[0094] (1.9)根据已经得到的X方向流速un+1、y方向流速，1和连续方程（1)，求得垂向 流速ω η+1，根据〇反坐标变换求得笛卡尔坐标系下垂向流速wn+1;
[0095] (1. 10)重复步骤（1.6)-(1. 9)，循环迭代计算，得到湖泊水域内不同位置和时间 的三维流场U、V、W及其静水位偏移位移η，并根据静水位偏移位移η进一步求得水深H。
[0096] (2)通过步骤（1)求得湖泊三维流场u、V、w和水深H等水动力参数的基础上，采 用分裂算法分别将水温控制方程和水质控制方程中的算子分为低频慢过程算子和高频快 过程算子两类，对其中的低频慢过程算子采用显式处理，对高频快过程算子采用隐式处理， 通过分裂算法分步离散并求解水温控制方程和水质控制方程，得到湖泊水域内不同位置和 时间的水温T和水质指标浓度C。
[0097] 具体而言，本步骤包括以下子步骤：
[0098] (2. 1)根据步骤（1. 5)类似的原理，将水温控制方程（5)和水质控制方程（6)的算 子分为低频慢过程算子和高频快过程算子两类，其中，低频慢过程算子包括对流项、水平扩 散项，高频快过程算子包括垂向扩散项；
[0099] (2. 2)对低频慢过程算子对流项和水平扩散项做显式处理，求得到中间温度Tl和 中间浓度Cl ;
[0100] (2. 3)根据求得的Tl、C1，对高频快过程算子垂向扩散项做隐式处理，求得下一个 时间节点的温度Tn+1、浓度C n+1，整理得到一个线性方程式，以温度T为例：
[0101]
[0102] M、S、P、R是包含已知变量Tl的参数，在各个网格写出式（12)，分别得到一个涉及 水温时空关系和水质时空关系的线性矩阵方程组，采用托马斯法求得τ η+1，并用同样的方法 求得cn+1;
[0103] (2. 4)重复步骤（2. 2)~（2. 3)，循环迭代计算，得到沙湖内不同位置和时间的水 温T和水质指标浓度C ;
[0104] (2. 5)根据求得的三维流场u、v、w、水深H、水温T和水质指标浓度C进行可视化 后处理，具体为：
[0105] (a)根据某时刻的湖区三维流场的计算结果生成流场矢量图（图6);
[0106] (b)对某个具体的网格，以时间为横坐标，以湖泊水环境某因子（流场、水温、水质 指标）的值为纵坐标，生成该因子的时间序列图（图7);
[0107] (C)根据某时刻的湖区垂向每层的水温（或水质）计算结果生成三维水温（或水 质）分布渲染图（图8)。
[0108] 模型验证
[0109] >计算值与实测值对比
[0110] 图9、图10