海洋数值模型的理论及应用

第一节海洋数值模型的发展概况观测，试验理论分析数值模型海洋的研究手段海洋数值模型发展的客观要求海水运动方程推导过程中的动力学假设，比如海水静力学假设，都是建立在大部分空间尺度大于10km，时间尺度大于惯性周期（2π/f）的情况下的海水运动方程不存在精确解：（1）偏微分方程本身是非线性的，（2）将一个实际问题具体化时所需要的环境场如海底地形、海岸线几何形状、表面驱动力等都无法给出解析函数表达式，只能在离散化的空间和时间区域取得因此，求解海水运动方程只能通过近似求解的办法来确定，而这些数值解又可能引入明显的近似误差海洋数值模型发展的历史海洋模型到按其水平网格的离散方式以及所使用的垂向坐标系的不同大致经历了如下几个发展阶段最早出现并且还在使用的海洋模型是Bryan等人开发的基于原始方程的低阶精度的有限差分模型，它在水深方向采用z坐标系目前常见的HOPS（HarvardOceanPredictionSystem），MOM（GFDLModularOceanModel），POP(ParallelOceanProgram)，NCOM（NCARCommunityOceanModel）模型在某种程度上都可以认为是该模型改进版海洋数值模型发展的历史上个世纪70年代，sigma坐标系开始应用于海洋模型在水深方向，比如目前被广泛使用的POM（PrincetonOceanModel）、ECOM（EstuarineCoastalandOceanModel）、ROMS（RegionalOceanModelingSystem）模型都属于这种类型的模型与传统的z坐标系相比，sigma坐标可以更好地贴合海底地形变化海洋数值模型发展的历史自上世纪90年代以来，非结构化网格技术开始在海洋模型中得到应用，相应的也出现了基于非结构化网格离散方式的有限元模型，如SEOM模型（SpectralFiniteElementOceanModel），和有限体积模型，如FVCOM虽然与传统的结构化网格相比,非结构化网格可以更好地拟合陆地边界，但是代码实现上的困难以及计算稳定性的问题使其迄今还没有得到非常广泛的应用新一代的海洋数值模型新一代的海洋模型广泛采用随地坐标系（terrain-followingcoordinates），进而促进了有关时间步长，对流项和压力梯度项等数值算法的改进进入新世纪以来，下一代的海洋数值动力模型正在紧锣密鼓的研制中，代表性的是TOMS(Terrain-followingOceanModelingSystem)，它融合了目前最先进的物理知识、数值方法和数据同化技术海洋模型的商业软件就海流的仿真建模研究而言，现阶段国外已形成了不少具有代表性的、具有强大前后处理功能及核心仿真建模技术的软件或系统（商业软件）如加拿大DalhousieUniversity研究的DALCOAST河口海岸预报系统；丹麦的DHIWater&Environment机构开发的的MIKE系列软件系统；荷兰的WLDelftHydraulics开发的的Delft－3D软件第二节基于POM模型的一般介绍水动力模型基本特点垂直方向的𝝈坐标变换水平网格采用正交曲线坐标和ArakawaC差分格式水平时间差分采用显格式，而垂直差分为隐格式自由表面可以模拟水位变化垂直和水平方向的混合扩散分别采用2.5阶的Mellor-Yamada湍流闭合模式和Smagorinski模式内外模态分别处理速度较慢的内重力波和速度较快的外重力波以提高整个模式计算效率包含了海水的热动力过程Sigma坐标变换：hs=-1z=0(a)s=0V(i,j+1)V(i,j)U(i+1,j)U(i,j)𝜼(i,j)(b)模型采用（a）垂向sigma坐标系和（b）水平ArakawaC网格zHhsh-=(1)控制方程0DUDVxyths=2UDUDUVDUfVDgDtxyxhs-20'''''MxoKgDDUdFxDxDssssss-=2VDUVDVDVfUDgDtxyyhs20'''''MyoKgDDVdFyDyDssssss-=连续方程海水运动方程(2)(3)(4)温度方程湍流动能方程HTKTDTUDTVDTTRFtxyDssss=-HSKSDSUDSVDSSFtxyDsss=盐度方程22222qKqDUqDVqDqqtxyDsss=3221222()()MHqoKUVgDqKFDBlsss-22222qKqlDUqlDVqlDqlqltxyDsss=322131()()MHloKUVgDqEEKWFDBsss-湍流混合长度方程(5)(6)(7)(8)其中，ω是基于sigma坐标系的垂向流速矢量，它与笛卡儿坐标系下的垂向流速W之间的关系可表示为dzzzzWUVdttxyDDDUVttxxyyhhhsss==='o=另外，方程中分别表示海水的实际密度，Boussinesq近似密度以及密度扰动其中，是实际水深DHh=(9)M—Y湍流闭合模型（profq）,,MHqKKK21(/)WElkL=在湍流动能及混合长度方程中分别表示湍流混合系数，热扩散系数和湍流动能的垂直扩散系数；墙近似函数，其中，κ=0.4是vonKarman常数湍流动能和混合长度方程组由下列等式关系闭合111()()LzHzh---=--21(/)WElL=0.2MMHHqKlqSKlqSKlq===其中，是稳定函数,MHSS(10)在2.5阶M—Y湍流闭合模型中，可表示为其中，。在不稳定层次条件下，的上限值为0.023，对应的值分别为2.0145和2.4401；在稳定层次条件下，它下限值为-0.28，对应的值分别为0.0470和0.0461。是2.5阶M-Y湍流闭合模式参数，由实验测得,MHSS0.42753.354134.67616.1270.494134.676hMhhHhGSGGSG-=--=-2201hglGqDs=10DshG10Ds1231,,,EEEB(11),MHSS,MHSSM—Y湍流模型的适用条件M-Y湍流模型已被广泛应用于浅海潮汐，风生表面混合层的模拟以及底边界层的研究；该模型在混合较弱的层化流体中对湍流混合系数计算效果不佳，应用于河口的可行性也有待于探讨；虽然该模型无论是在物理上还是在数学上都存在着明显的缺陷，但是目前依然得到广泛应用。只有深入了解海水的微细结构以及海洋湍流结构和性质后，才有可能构建出更为完善的湍流闭合模型海水运动方程的水平扩散项•其中水平切应力•是水平涡粘性系数，由Smagorinsky模式计算()()()()xxxxyyxyyyFHHxyFHHxy==2()2xxMxyyxMyyMUAxUVAyxVAy====MA(12)(13)Smagorinsky水平扩散模式•其中参数C为无量纲值，其取值为0.1或0.2，如果计算网格划分得足够细的话，C可以取值为0•由Smagorinsky模式计算的水平涡粘性系数随着网格精度的改善和速度梯度的减小而减小1/222211=22MuuvvACxyxxyy(14)温度盐度的水平扩散项•其中•是水平热扩散系数，一般地，是一个小值，取0.1或0.2，甚至在某些情况下可以取为0HA22()(),,,,xyFHqHqTSqqlxy==22xHyHqAxqAy==HMAA(15)(16)垂向边界条件1/2221010101022212/3,(,)(0)0,(0),(0)(0),(0)=(0),0MasHKuvCUVUVDKTSwwsDqqBussss=-==-1.海面0s(17)2.海底1s-1/222222/321,(,)(1)0,0(1),(1)=(1),0MbbbbbHKuvCuvuvDKTSDqqBussss=-==---其中，分别是湍流闭合模型参数和摩擦速度；1,Bu,sbCC分别是海面和海底的摩擦系数；(18)221=max,0.0025ln1+/bkboCHzs-模型的外模态•控制近岸环流的海水运动方程包含了运动速度较快的外重力波和速度较慢的内重力波，因此进行模式分裂可以大大提高模型的计算效率•模型的外模态就是计算自由水位和垂向平均的流速变化，需要较小的时间步长•内模态主要模拟流速矢量、温度、盐度等三维结构变化，对时间步长的要求相对宽松•采用模式分裂技巧可以有效地减少因外模态计算所消耗的计算时间，从而提高计算效率垂向积分的海水运动方程0DUDVxyth=2(0)+(-1)xUDUDUVDFfVDgDwuwutxyxh--=--1''+ooxogDDGDddxxsssss--2(0)+(-1)yVDUVDVDFfUDgDwvwvtxyyh-=--1''+ooyogDDGDddyysssss---1-1,ooUUdVVdss==其中：连续方程海水运动方程(19)(20)(21)水平扩散项=2++=++2xMMyMMUUVFHAHAxxyyxUVVFHAHAxyxyy水平数值耗散项2222xxxyyyUDUVDUDUVDGFFxyxyUVDVDUVDVDGFFxyxy=---=---(22)(23)模型代码中的一些符号•im,jm,kb,imm1,jmm1,kbm1,mode,isplit•dte,dti,days,umol•z,zz,dz,dzz•aam2d,art,aru,arv,dum,dvm,fsm•dx,dy,h,el,d,ua,va,ut,vt,et,cor•swrad•wusurf,wvsurf,wubot,wvbot,wtsurf,wssurf,rad•u,v,ω,t,s,rho,l,q2•km,kh,kq,aam,aah•rmean,tclim,sclim数值积分流程图START9000IINT=1,IENDPrintSTOPSetParmeters,InitialValuesSetParametersInitialValuesBAROPG8000AdjustintegralofU,VtomatchUT,VTVERTVLBCOND(5)ADVQ(Q2)ADVQ(Q2L)PROFQBCOND(6)ADVT(T)ADVT(S)PROFT(T)PROFT(S)BCOND(4)ADVUADVVPR