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

第一节海洋数值模型的发展概况观测，试验理论分析数值模型海洋的研究手段海洋数值模型的理论及应用海洋数值模型发展的客观要求海水运动方程推导过程中的动力学假设，比如海水静力学假设，都是建立在大部分空间尺度大于10km，时间尺度大于惯性周期（2π/f）的情况下的海水运动方程不存在精确解：（1）偏微分方程本身是非线性的，（2）将一个实际问题具体化时所需要的环境场如海底地形、海岸线几何形状、表面驱动力等都无法给出解析函数表达式，只能在离散化的空间和时间区域取得因此，求解海水运动方程只能通过近似求解的办法来确定，而这些数值解又可能引入明显的近似误差海洋数值模型的理论及应用海洋数值模型发展的历史海洋模型到按其水平网格的离散方式以及所使用的垂向坐标系的不同大致经历了如下几个发展阶段最早出现并且还在使用的海洋模型是Bryan等人开发的基于原始方程的低阶精度的有限差分模型，它在水深方向采用z坐标系目前常见的HOPS（HarvardOceanPredictionSystem），MOM（

GFDL

ModularOceanModel），POP

(ParallelOceanProgram)，NCOM（NCARCommunityOceanModel）模型在某种程度上都可以认为是该模型改进版海洋数值模型的理论及应用海洋数值模型发展的历史上个世纪70年代，sigma坐标系开始应用于海洋模型在水深方向，比如目前被广泛使用的POM（PrincetonOceanModel）、ECOM（EstuarineCoastalandOceanModel）、ROMS（RegionalOceanModelingSystem）模型都属于这种类型的模型与传统的

z

坐标系相比，sigma坐标可以更好地贴合海底地形变化海洋数值模型的理论及应用海洋数值模型发展的历史自上世纪90年代以来，非结构化网格技术开始在海洋模型中得到应用，相应的也出现了基于非结构化网格离散方式的有限元模型，如SEOM模型（SpectralFiniteElement

Ocean

Model），和有限体积模型，如FVCOM虽然与传统的结构化网格相比,非结构化网格可以更好地拟合陆地边界，但是代码实现上的困难以及计算稳定性的问题使其迄今还没有得到非常广泛的应用海洋数值模型的理论及应用新一代的海洋数值模型新一代的海洋模型广泛采用随地坐标系（terrain-followingcoordinates），进而促进了有关时间步长，对流项和压力梯度项等数值算法的改进进入新世纪以来，下一代的海洋数值动力模型正在紧锣密鼓的研制中，代表性的是TOMS(Terrain-followingOceanModelingSystem)，它融合了目前最先进的物理知识、数值方法和数据同化技术

/WWWPUBLIC/htdocs.pom/TOMS.htm

海洋数值模型的理论及应用海洋模型的商业软件就海流的仿真建模研究而言，现阶段国外已形成了不少具有代表性的、具有强大前后处理功能及核心仿真建模技术的软件或系统（商业软件）如加拿大DalhousieUniversity研究的DALCOAST河口海岸预报系统；丹麦的DHIWater&Environment机构开发的的MIKE系列软件系统；荷兰的WLDelftHydraulics开发的的Delft－3D软件海洋数值模型的理论及应用第二节

基于POM模型的一般介绍海洋数值模型的理论及应用水动力模型基本特点垂直方向的

𝝈坐标变换水平网格采用正交曲线坐标和ArakawaC差分格式水平时间差分采用显格式，而垂直差分为隐格式自由表面可以模拟水位变化垂直和水平方向的混合扩散分别采用2.5阶的Mellor-Yamada湍流闭合模式和Smagorinski模式内外模态分别处理速度较慢的内重力波和速度较快的外重力波以提高整个模式计算效率包含了海水的热动力过程海洋数值模型的理论及应用Sigma坐标变换：hs=-1z=0z=H(x,y)

s=0

V(i,j+1)

V(i,j)

U(i+1,j)

U(i,j)

𝜼

(i,j)模型采用（a）垂向sigma坐标系和（b）水平ArakawaC网格(1)海洋数值模型的理论及应用控制方程连续方程海水运动方程(2)(3)(4)海洋数值模型的理论及应用温度方程湍流动能方程盐度方程湍流混合长度方程(5)(6)(7)(8)海洋数值模型的理论及应用其中，ω是基于sigma坐标系的垂向流速矢量，它与笛卡儿坐标系下的垂向流速W之间的关系可表示为另外，方程中分别表示海水的实际密度，Boussinesq近似密度以及密度扰动其中，是实际水深(9)海洋数值模型的理论及应用M—Y湍流闭合模型（profq）在湍流动能及混合长度方程中分别表示湍流混合系数，热扩散系数和湍流动能的垂直扩散系数；墙近似函数，其中，κ=0.4是von

Karman常数湍流动能和混合长度方程组由下列等式关系闭合其中，是稳定函数(10)海洋数值模型的理论及应用在2.5阶M—Y湍流闭合模型中，可表示为其中，。在不稳定层次条件下，的上限值为0.023，对应的值分别为2.0145和2.4401；在稳定层次条件下，它下限值为-0.28，对应的值分别为0.0470和0.0461。

是2.5阶M-Y湍流闭合模式参数，由实验测得(11)海洋数值模型的理论及应用M—Y湍流模型的适用条件M-Y湍流模型已被广泛应用于浅海潮汐，风生表面混合层的模拟以及底边界层的研究；该模型在混合较弱的层化流体中对湍流混合系数计算效果不佳，应用于河口的可行性也有待于探讨；虽然该模型无论是在物理上还是在数学上都存在着明显的缺陷，但是目前依然得到广泛应用。只有深入了解海水的微细结构以及海洋湍流结构和性质后，才有可能构建出更为完善的湍流闭合模型海洋数值模型的理论及应用海水运动方程的水平扩散项其中水平切应力是水平涡粘性系数，由Smagorinsky模式计算(12)(13)海洋数值模型的理论及应用Smagorinsky水平扩散模式其中参数C

为无量纲值，其取值为0.1或0.2，如果计算网格划分得足够细的话，C可以取值为0

由Smagorinsky模式计算的水平涡粘性系数随着网格精度的改善和速度梯度的减小而减小(14)海洋数值模型的理论及应用温度盐度的水平扩散项其中是水平热扩散系数，一般地，是一个小值，取0.1或0.2，甚至在某些情况下可以取为0(15)(16)海洋数值模型的理论及应用垂向边界条件1.海面(17)海洋数值模型的理论及应用2.海底其中，分别是湍流闭合模型参数和摩擦速度；分别是海面和海底的摩擦系数；(18)海洋数值模型的理论及应用模型的外模态控制近岸环流的海水运动方程包含了运动速度较快的外重力波和速度较慢的内重力波，因此进行模式分裂可以大大提高模型的计算效率模型的外模态就是计算自由水位和垂向平均的流速变化，需要较小的时间步长内模态主要模拟流速矢量、温度、盐度等三维结构变化，对时间步长的要求相对宽松采用模式分裂技巧可以有效地减少因外模态计算所消耗的计算时间，从而提高计算效率海洋数值模型的理论及应用垂向积分的海水运动方程其中：连续方程海水运动方程(19)(20)(21)海洋数值模型的理论及应用水平扩散项水平数值耗散项(22)(23)海洋数值模型的理论及应用模型代码中的一些符号im,jm,kb,imm1,jmm1,kbm1,mode,isplitdte,dti,days,umolz,zz,dz,dzzaam2d,art,aru,arv,dum,dvm,fsmdx,dy,h,el,d,ua,va,ut,vt,et,corswradwusurf,wvsurf,wubot,wvbot,wtsurf,wssurf,radu,v,ω,t,s,rho,l,q2km,kh,kq,aam,aahrmean,tclim,sclim海洋数值模型的理论及应用数值积分

流

程

图海洋数值模型的理论及应用外模态和内模态的相互作用海洋数值模型的理论及应用内模态的计算方法三维变量的计算（以T为例）可以分解为垂直扩散项的计算和水平对流扩散项的计算，以温度方程为例方程写为求解过程分为两步，第一步计算水平对流扩散项，采用中心差分格式，由advt子程序实现(24)(25)水平对流扩散项垂直扩散项海洋数值模型的理论及应用第二步计算垂向扩散项，数值方法采用“蛙跳”格式，由proft子程序实现由于“蛙跳”格式在奇数时间步时的解与在偶数时间步时的解会发生偏离，因此每一时间步的计算结束后还需对上述解进行平滑处理，即α为一个小值，一般取0.05。处理后，(26)(27)海洋数值模型的理论及应用计算网格的安排二维外模态网格分布三维内模态网格分布海洋数值模型的理论及应用水平对流项的计算■虽然模型采用有限差分计算格式，但是对于每个网格对流项的计算都是按照有限容积的方法进行处理，即温度对流过程可表示为■速度的对流过程与温度相似，可写为其中，是由z坐标转化为σ坐标后产生的弯曲项(28)(29)海洋数值模型的理论及应用垂直扩散项的计算■垂直扩散项的计算公式（第k层，1<k<kb-1）■可以改写为(30)(31)海洋数值模型的理论及应用■展开后合并同类项得到其中(32)(33)海洋数值模型的理论及应用■现在假定温度解的形式为■由(34)得到的代入到(32)就可以得到其中的值由(33)求得，的值由上一层的值确定(34)(35)海洋数值模型的理论及应用■当k=1时，即表层海水，海水温度主要取决于海表面温度通量。公式（34）的解可以近似写为■这样，表层海水的温度只能根据式（31）求解，即■

上式可以进一步表示为(36)(37)海洋数值模型的理论及应用■式（37）还可以进一步写为■再与温度的通解（34）比较■

可以得到(34)(38)(39)海洋数值模型的理论及应用■当短波辐射通量的计算公式■

r，ad1，ad2均为光辐射常数，根据不同水质取值如下(40)ntp

12345JerlovtypeIIaIbIIIIIr

0.580.620.670.700.78ad1(m)0.350.60

1.01.51.4ad2(m)

23.020.017.014.07.9来源：Jerlov，1976；PaulsonandSimpson，1977海洋数值模型的理论及应用■当k=kb-1时，即底层海水，假定海底的热通量为0，根据前面的推导方式同样可以得到底层海水温度的计算公式■

对于盐度方程垂向扩散项的计算与温度方程相同，只是不考虑太阳短波辐射这一项(41)海洋数值模型的理论及应用时间步长的CFL限制条件●由于模型水平方向采用显格式，因此时间步长的选取必须要满足CFL稳定条件。●对于外模态，时间步长的限制条件为其中，可能预见到的最大流速●对于内模态，时间步长的限制条件较外模态的情形宽松很多，主要是速度较快的外重力波已经在外模态中考虑了。一般取值30~50即可

(42)海洋数值模型的理论及应用时间步长的其它限制条件●对于动量或标量还有其它一些时间限制条件其中或●以及

其中分别为地球角速率和地理纬度(43)(44)海洋数值模型的理论及应用侧开边界条件（bcond）陆地及岸线是由dum、dvm、fsm控制的，在陆地上这些变量的值设为0，有水的地方设为1子程序bcond（idx），idx=1对应的水位边界条件；idx=2对应垂向平均流速边界条件；idx=3对应三维水平流速边界条件；idx=4对应温盐边界条件；idx=5对应垂向流速边界条件；idx=6对应湍流动能和湍流混合长度边界条件模型对开边界条件的要求很高，而开边界条件本身又具有很大的不确定性，因此有必要对模型的外模态和内模态的开边界分别进行处理海洋数值模型的理论及应用数值积分

流

程

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

EASTuaf(im,j)=2*bc(j)/(h(im,j)+elf(im,j)+h(imm1,j+elf(imm1,j))elf(im,j)=elf(imm1,j)vaf(im,j)=setWESTuaf(2,j)=2*bc(j)/(h(1,j)+elf(1,j)+h(2,j)+elf(2,j))elf(1,j)=elf(2,j)vaf(1,j)=setNORTHvaf(i,jm)=2*bc(i)/(h(i,jm)+elf(i,jm)+h(i,jmm1)+elf(i,jmm1))elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTHvaf(i,2)=2*bc(i)/(h(i,1)+elf(i,1)+h(i,2)+elf(i,2))elf(i,1)=elf(i,2)uaf(i,1)=setElevationconditionh

=BC

EASTelf(imm1,j)=bc(j)elf(im,j)=elf(imm1,j)cosmeticuaf(im,j)=uaf(imm1,j)vaf(im,j)=setWESTelf(2,j)=bc(j)uaf(2,j)=uaf(3,j)vaf(1,j)=setNORTHelf(i,jmm1)=bc(i)elf(i,jm)=elf(i,jmm1)cosmeticvaf(i,jm)=vaf(i,jmm1)uaf(i,jm)=setSOUTHelf(i,2)=bc(i)vaf(i,2)=vaf(i,3)uaf(i,1)=set

外模态

开

边

界

条

件

一海洋数值模型的理论及应用FormulaBoundaryCodeRadiationEASTuaf(im,j)=sqrt(grav/h(imm1,j))*el(imm1,j)+bc(j)elf(im,j)=elf(imm1,j)vaf(im,j)=setWESTuaf(2,j)=-sqrt(grav/h(2,j))*el(2,j)+bc(j)elf(1,j)=elf(2,j)vaf(1,j)=setNORTHvaf(i,jm)=sqrt(grav/h(i,jmm1))*el(i,jmm1)+bc(i)elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTHvaf(i,2)=-sqrt(grav/h(i,2))*el(i,2)+bc(i)elf(i,1)=elf(i,2)uaf(i,1)=setRadiation

EASTgae=dte*sqrt(grav*h(im,j))/dx(im,j)uaf(im,j)=gae*ua(imm1,j)+(1.-gae)*ua(im,j)elf(im,j)=elf(imm1,j)vaf(im,j)=setWESTgae=dte*sqrt(grav*h(2,j))/dx(2,j)uaf(2,j)=gae*ua(3,j)+(1.-gae)*ua(2,j)elf(1,j)=elf(2,j)vaf(1,j)=setNORTHgae=dte*sqrt(grav*h(i,jm))/dy(i,jm)vaf(i,jm)=gae*va(i,jmm1)+(1.-gae)*va(i,jm)elf(i,jm)=elf(i,jmm1)uaf(i,jm)=setSOUTHgae=dte*sqrt(grav*h(i,2))/dy(i,2)vaf(i,2)=gae*va(i,3)+(1.-gae)*va(i,2)elf(i,1)=elf(i,2)uaf(i,1)=set

外模态

开

边

界

条

件

二流速水位流速海洋数值模型的理论及应用FormulaBoundaryCodeRadiationEASTgae=dte*sqrt(grav*h(imm1,j))/dx(imm1,j)elf(imm1,j)=gae*el(imm2,j)+(1.-gae)*el(imm1,j)elf(im,j)=elf(imm1,j)uaf(im,j)=uaf(imm1,j)vaf(im,j)=setWESTgae=dte*sqrt(grav*h(2,j))/dx(1,j)elf(2,j)=gae*el(3,j)+(1.-gae)*el(2,j)uaf(2,j)=uaf(3,j)vaf(2,jm)=setNORTHgae=dte*sqrt(grav*h(i,jmm1))/dy(i,jmm1)elf(i,jmm1)=gae*el(i,jmm2)+(1.-gae)*el(i,jmm1)elf(i,jm)=elf(i,jmm1)vaf(i,jm)=vaf(i,jmm1)uaf(i,jm)=setSOUTHgae=dte*sqrt(grav*h(i,2))/dy(i,2)elf(i,2)=gae*el(i,3)+(1.-gae)*el(i,2)vaf(i,2)=vaf(i,3)，uaf(i,1)=setCyclicEAST(I=IM)elf(im,j)=elf(3,j)uaf(im,j)=uaf(3,j)，vaf(im,j)=vaf(3,j)WEST(I=1)elf(1,j)=elf(imm2,j)，elf(2,j)=elf(imm1,j)uaf(2,j)=uaf(imm1,j)，vaf(2,j)=vaf(imm1,j)NORTH(J=JM)elf(i,jm)=elf(i,3)uaf(i,jm)=uaf(i,3)vaf(i,jm)=vaf(i,3)SOUTH(J=1)elf(i,1)=elf(i,jmm2)，elf(i,2)=elf(i,jmm1)uaf(i,2)=uaf(i,jmm1)vaf(i,2)=vaf(i,jmm1)

外模态

开

边

界

条

件

三水位海洋数值模型的理论及应用FormulaBoundaryCodeInflowconditionEASTuf(im,j,k)=bc(j,k)vf(im,j,k)=setU=BC

WESTuf(2,j,k)=bc(j,k)vf(1,j,k)=set

NORTHvf(i,jm,k)=bc(i,k)uf(i,jm,k)=setSOUTHvf(i,2,k)=bc(i,k)uf(i,1,k)=setRadiation:EASTgai=sqrt(h(im,j)/hmax)uf(im,j,k)=gai*u(imm1,j,k)+(1.-gai)*u(im,j,k)vf(im,j,k)=setWESTgai=sqrt(h(2,j)/hmax)uf(2,j,k)=gai*u(3,j,k)+(1.-gai)*u(2,j,k)vf(1,j,k)=setNORTHgai=sqrt(h(i,jm)/hmax)vf(i,jm,k)=gai*v(i,jmm1,k)+(1.-gai)*v(i,jm,k)uf(i,jm,k)=setSOUTHgai=sqrt(h(i,2)/hmax)vf(i,2,k)=gai*v(i,3,k)+(1.-gai)*v(i,2,k)uf(i,1,k)=set

内模态开边界条件（一）海洋数值模型的理论及应用FormulaBoundaryCodeUpstreamadvectiononTorSEASTuf(im,j,k)=t(im,j,k)-dti/(dx(im,j)+dx(imm1,j))*((u(im,j,k)+abs(u(im,j,k)))*(t(im,j,k)-t(imm1,j,k))+(u(im,j,k)-abs(u(im,j,k)))*(tbe(j,k)-t(im,j,k)))WESTuf(1,j,k)=t(1,j,k)-dti/(dx(1,j)+dx(2,j))*((u(1,j,k)+abs(u(1,j,k)))*(t(1,j,k)-tbw(j,k))+(u(1,j,k)-abs(u(1,j,k)))*(t(2,j,k)-t(1,j,k)))

NORTHuf(i,jm,k)=t(i,jm,k)-dti/(dy(i,jm)+dy(i,jmm1))*((v(i,jm,k)+abs(v(i,jm,k)