Challenges in ocean model development

1. Challenges in ocean model development Mats Bentsen Nansen Environmental and Remote Sensing Center Bjerknes Centre for Climate Research Nansen-Zhu International Research Centre

2. Outline • Model purpose • Vertical coordinate • Horizontal grid • Time stepping • Transport • Status and summary

3. Model purpose • Ocean general circulation model • The model should perform well on the global scale • as well as on regional scale • It should be suited to be a component of an • Earth System Model

4. Atmosphere Atmosphere Atmosphere Atmosphere Atmosphere Atmosphere Land surface Land surface Land surface Land surface Land surface Ocean & sea-ice Ocean & sea-ice Ocean & sea-ice Ocean & sea-ice Sulphate aerosol Sulphate aerosol Sulphate aerosol Non-sulphate aerosol Non-sulphate aerosol Carbon cycle Carbon cycle Atmospheric chemistry Sulphur cycle model Non-sulphate aerosols Ocean & sea-ice model Land carbon cycle model Carbon cycle model Ocean carbon cycle model Atmospheric chemistry Atmospheric chemistry EARTH SYSTEM MODEL DEVELOPMENT 1975 1985 1992 1997 2002 2005 Prognostic model development Strengthening colours denote improvements in models Off-line model development Strengthening colours denote improvements in models

5. What vertical coordinate to use?

6. Hybrid models • HYCOM: HYbrid Coordinate Ocean Model, University of Miami • HYPOP: Hybrid version of POP (Parallel Ocean Program), Los Alamos

7. Horizontal grids Structured grids Method: finite difference Unstructured grids Methods: finite element, spectral element

8. Staggered grids Arakawa A Arakawa B Arakawa C h u h h,u,v v u,v The properties of the grids can be studied with the linearized shallow- water equations

9. Assuming wave solutions give the dispersion relation where we have introduced the radius of deformation Given grid resolution d, we will study the σ/f as a function of kd and ld for two cases of λ/d.

10. Well resolved radius of deformation (λ/d = 2) Continous A-grid B-grid C-grid

11. Not resolved radius of deformation (λ/d = 0.1) Continous A-grid B-grid C-grid

12. Why does the various grids perform the way they do? Arakawa A: averaging in the evaluation of the height gradients h,u,v h,u,v h,u,v Arakawa B: averaging in the evaluation of the height gradients Arakawa C: averaging in the evaluation of the Coriolis terms h h u u h h v v u,v u,v h h u u h h u,v u,v v v

13. Recent methods to overcome the problems of traditional grids • C-D grid by Adcroft et al. (1999). • Shifted approximations by Nechaev and Yaremchuk (2004). • Reversibly staggered grids by McGregor (2005).

14. C-D grid C-grid D-grid C-D-grid • Solves for velocities for both C- and D-grid • No need for averaging in the evaluation of Coriolis terms • A computational mode is introduced because twice as many velocity • variables is used compared to the C-grid • The computational mode is damped by using an implicit backward • scheme for the Coriolis terms.

15. Shifted approximations u u h h v v u u h h v v • A C-grid is used where the Coriolis term is approximated at 4 • shifted locations • The solution is obtained by averaging the shifted approximation • The scheme introduces damping of waves with high wavenumbers

16. Reversibly staggered grid (R-grid) h,u,v • Variables are stored as in the A-grid • An interpolation is used that preserves amplitude, but have a phase • error for short waves • Velocities are interpolated to C-grid positions for the computation of • the divergence • Height gradients are evaluated at C-grid locations and interpolated to • A-grid locations

17. Not resolved radius of deformation (λ/d = 0.1) R-grid “left” R-grid “right”

18. Boundary pressure adjustment experiment by Milliff and McWilliams (1994). The initial condition is a nearly balanced vortex that evolves in a non-linear manner and interacts with the boundary. The dynamic pressure is shown. C-grid solution after 200 days R-grid solution after 200 days

19. Time stepping We will consider the 1D hyperbolic system of equations which is an simple analog to the barotropic mode in the absence of Coriolis force and topography. Assuming a wave solution, we can express the evolution of Fourier components for wavenumber k as In the case of the Leapfrog algorithm, these equations can be discretized as where

20. The Leapfrog algorithm has no amplitude error but rather big • phase error • It has an computational mode that need to be filtered • Application of the traditional Asselin time filter increases the phase • error, causes considerable damping of well resolved waves, and • reduces the maximum allowable time step Roots of the characteristic equation of Leapfrog time-stepping

21. Shchepetkin and McWilliams (2005) studied a range of generalized forward-backward algorithms. The forward-backward algorithm can be written It has a stability limit of α = 2, no amplitude error, but quite large phase error. Roots of the characteristic equation of forward-backward time-stepping

22. A generalized forward-backward algorithm has a predictor step followed by a corrector step The particular weights gives a 3rd order accurate algorithm with stability limit α ≈ 2.14, damping of short waves and a small phase error.

23. Roots of the characteristic equation of the generalized forward-backward time-stepping

24. Transport Proper treatment of transport processes is very important in many geophysical applications. We are often faced with a continuity equation together with a tracer equation • Desired properties of transport algorithms: • Accuracy in terms of amplitude, shape, and phase • Conservation • Compability • Efficiency

25. 2D advection test (split cylinder) Solid body rotation as described in Zalesak, J. Comput. Phys, 31, 1979. To complete one revolution, 628 time steps are needed. Initial condition

26. MPDATA [0…1] MPDATA [1…2]

27. FCT2 FCT4

28. Weighted essentially non-oscillatory (WENO) scheme

29. WENO3

30. WENO3 WENO5

31. Incremental remapping (Dukowicz and Baumgardner, 2000)

32. With remapping there is a small incremental cost of adding more tracers, since many calculations used to update the first tracer can be reused • Incremental remapping has the properties: • Accuracy in terms of amplitude, shape, and phase • Conservation • Compability • Efficiency (for many tracers)

33. REMAP

34. Transport of temperature and salinity in an isopycnic layer • In the continuous case: • Isopycnic advection does not lead to density change • Isopycnic diffusion increases density (cabbeling)

35. Status and summary • Our development concentrates on the use of isopycnic vertical • coordinate, with possible extension to a hybrid coordinate. • We have not decided on which horizontal grid to use, and active • research is going in to this. At the moment, the R-grid looks • promising but needs testing in more realistic situations. • The generalized forward-backward algorithm shown, or a variant • of it will most likely be used. • We will most likely use incremental remapping as the transport • algorithm. It is a method suited for B-grid models but we have • adapted it to be used with C-grid