Numerical Aspects of Marine Biogeochemical Modelling

1. Numerical Aspects of Marine Biogeochemical Modelling

2. Outline Continuum hypothesis and balance equations: When is the continuous approximation of a discrete world appropriate? The tracer transport equation. Discretisation of PDEs: from mathematical equations to discrete numerics. Concepts of convergence, consistency and discretisation error, mass conservation, positivity. Time integration methods: Single-step methods Mutlistep-methods Implicit vs. explicit schemes The advection-diffusion operator Advection schemes Diffusion schemes Coupling methods for biogeochemical dynamics and the advection-diffusion process Why splitting? Splitting Methods (conceptually) Splitting Methods mathematically � splitting error The role of Scales Introduction to the model used in the exercises

3. Marine Biogeochemical Models The type of models considered in this lecture, as the vast majority of biogeochemical models, quantify concentration changes in given locations. are Eulerian (y quantity at fixed positions in space) as oposed to Lagrangian (y attribute of a moving particle). are continuous.

4. Continuum Hypothesis Any quantity in the model is a continuous function of space and time, i.e. to any location in time and space can be attributed a value of this quantity and the changes are continuous. To apply the continuum hypothesis in a world that is essentially discrete, we must look at reality from a distance such that discontinuities smooth out.

5. Scales and the validity of the continuum hypothesis Microscopic range: measurements show non-continuous fluctuations. Mesoscopic range: measurements show constant values. Macroscopic range: measurements show continuous variations.

6. Scales and the validity of the continuum hypothesis

7. The Balance Equation A quantity in a volume may change through: Advection Efflux (in relation to surrounding) Supply (independent of surrounding This concept is generally applicable: Mass balance Conservation of momentum Conservation of energy Navier-Stokes

8. The Balance Equation

9. The tracer transport equation

10. Discretisation Necessary for numerical treatment. Aim: replacing the continuous space-time domain with a discrete grid domain, in order to compute integrals and derivatives with algebraic expressions.

11. How to translate the equations into discrete form? Discretisation through Finite differences Finite volume Finite elements

12. Finite Difference Method A derivative is defined as Replace derivative equations with corresponding difference quotient

13. Finite Volume Method The domain is discretised by small cells in which all quantities are supposed to be distributed homogenously and fluxes are defined on the cell faces. => based on the finite volume formulation of the balance equation.

14. Important numerical concepts � Discretisation Error Errors induced by the approximation of a continuous model by a discrete system. Other errors e.g. Model approximation errors (conceptual, included in the analytical solution). Rounding errors (numerical, 0->0.000000001).

15. For decreasing step-sizes the equations of the discrete system approach the original pde. Based on local error. => Evaluated on single time step Important numerical concepts - Consistency

16. Important numerical concepts � Stability and Stiffness A system is stable when its solution is not sensitive to small perturbations in the initial or boundary conditions. A discretisation is stable if it doesn�t amplify small perturbations such as rounding errors. A pde-system is called stiff, when it involves modes on very different scales that challenge the stability of the discrete solution (while the true solution is barely unaffected).

17. Important numerical concepts - Convergence For decreasing step sizes the discrete system solution approaches an asymptotic value, which is the true solution of the system. Related to global error.

18. Discrete Time integration Single-step schemes Euler Runge-Kutta Multistep schemes Adams Leap-frog

19. The Taylor Theorem Given a function and its derivatives are finite and continuous, the value at each point t+Dt is given by: Cut-off error is increasingly small for Dt < 1. Extremely useful for determining discretisation errors. Designing discretisation schemes.

20. Single Step Methods Based only on current time-step. Higher order schemes require multiple evaluation of RHS- function. Memory friendly. Relatively simple to implement.

21. Euler Scheme: Linear extrapolation Simple, Fast. 1st order.

22. Runge-Kutta Schemes Idea: evaluating F on multiple intermediate positions of Dt to obtain a higher order approximation. Ex.: Runge-Kutta of 2nd order (Heun-method)

23. Runge-Kutta 4th order Precise to the 4th order Robust Widely used in commercial ode solvers

24. Multi Step Schemes Involve previous time steps. Don�t require additional evaluations. Require more memory. Are slightly more tedious to implement, especially when splitting methods are involved.

25. Adams Schemes Idea: approximation of F by polynomial extrapolation of previous timesteps. Ex. Adams-Bashforth 2nd order (linear extrapolation): Weakly unstable.

26. Leap-frog Idea: interpolating linearly from mean previous change.

27. Asselin-Filter Numerical filter for instability waves, mostly used to smooth leap-frog solutions.

28. Mass conservation In a closed boundary system the total mass has to be conserved.

29. Example: Simple NPD model

30. A stiff system

31. Explicit vs. implicit schemes Schemes considered so far are explicit, i.e. contain only known values on the RHS:

32. Predictor-Corrector Method Simple way to approximate a implicit solution: Solve the system explicitly. (predictor step) Insert the computed �future solution� in an implicit solver Ex.: Runge-Kutta method from before.

33. Should I decrease my time step unlimited if I am not restricted by computational resources? No!

34. Advection schemes Upwind schemes Central schemes MPDA TVD schemes, e.g. MUSCL Piecewise parabolic schemes

35. Upwind scheme The flux on the cell faces is approximated by the upstream concentration:

36. Central scheme The flux on the cell faces is approximated by the mean of the adjacent cells:

37. Concept of numerical diffusion in advection schemes The discretisation error may add an artificial numerical diffusion term to the original equations. Smoothes strong gradients. Enhances stability, but gives wrong results.

38. Concept of dispersion in advection schemes Introduction of spurious oscillations to the solution, particularly at strong curvatures. Scheme introduces numerical modes that travel at frequency dependent speed. Can lead to instabilities of the solution. It can be shown that a non dispersive-scheme can be no higher than first-order.

39. Concept of positivity in advection schemes The advected field shall not become negative at any point in space or time! Achieved by the construction of schemes that maintain the advected field in the current range of total variation.

40. Upwind scheme Non dispersive Diffusive Positive Monotonicity preserving 1st order

41. Central scheme 2nd order Non-diffusive Dispersive Unconditionally unstable

42. Numerical diffusion in upwind scheme Applying Taylor series in time and space to the upwind scheme leads to:

43. Upwind vs. Central Scheme

44. Multidimensional Positive Definite Advection Scheme Family of upwind schemes Numerical diffusion is corrected through antidiffusive velocities. Correction only applied when positivity is guaranteed.

45. Total Variance Diminishing Schemes Blending of upwind and central scheme. Blending regulated by flux limiter that controls the total variance of a scheme. Monotonicity conserving. Posivtive

46. Monotonic Upstream Centred Scheme for Conservation Laws (MUSCL) Slope interpolation of cell centre values. Slope is limited by appropriate TVD criterion. Positive Monotonicity conserving. High order when unlimited.

47. Piecewise parallel scheme Piecewise parallel interpolation of cell centre points. Flux limiter determines curvature. 3rd order. Positive, monotonicity conserving.

48. Diffusion Usually central 2nd order schemes are applied. Vertical diffusion involves strong mixed layer turbulence. This in combination with the small vertical dimension of the grid cells favours implicit solvers for the vertical mode.

49. Coupled Systems Systems composed of several sub-modules. Sub-Modules can be integrated as a whole or split. Why splitting? Why not splitting? Methods: Operator Splitting Source Splitting

50. Why splitting a coupled systems? Different processes may have different scales and therefore different discretisation errors and stability requirements. E.g.: Advection + hor. diffusion, then vertical diffusion. Other ex.: Atmospheric chemistry. Memory Combination of spatial and time derivates favours certain combined choices, e.g.: Euler Forward + upwind family, leap-frog + centred schemes with lagged diffusion�

51. Operator Splitting Method Two processes applied successively. Straight forward to implement. Process integration is decoupled => spurious transition. First order.

52. Operator Splitting Method

53. Source Splitting Method The supposingly weaker varying processes is estimated piecewise-linear on the global time step. Consistent. First order.

54. Source Splitting Method

55. Why not splitting? The uncoupling of the two processes involved causes a splitting error as the variation of one process is not considered over the integration of the other. Taylor expansions show that for both methods the error is of 1st order.

56. Problems Which is the dominant scale? Which is the dominant error?

57. Implementation Relatively easy and numerically cost-efficient Largely independent of sub-process integrations Not always, e.g.: If the weaker varying process is estimated with a leap-frog scheme, the estimation has to be accordingly:

58. Example: 1D fully coupled biogeochemical model Bio geochemical model: BFM Hydrodynamical model: POM 1D set-up in Northern Adriatic Investigation: Does the biogeochemistry need a higher resolution than the physics? Does the time integration scheme used in split mode have an influence on the solution? Which error dominates the solution? Which scale dominates the solution?

59. Does the biogeochemistry need a higher resolution than the physics? Errors at two time steps differing by more than one order of magnitude against a reference solution. No significant difference in solution => No high resolution of biogeochemistry needed.

60. Does the time integration scheme used in split mode have an influence on the solution? Various integration schemes estimating the physical rate in source splitting. No significant difference in solution => Integration scheme has no influence.

61. Which error dominates the solution? Source splitting estimating the physics against full integration. Splitting provokes significant errors that dominate. Operator Splitting produces similar results.

62. Which scale dominates? Scale of mixed layer turbulence processes dominates over biogeochemical time scale.

63. NPD2 � A simple NPD model with two vertical boxes Two part exercise: Spin-up the two boxes to steady state, no diffusion. Fully coupled dynamics with turbulence cycle.

64. NPD2 � A simple NPD model with two vertical boxes Biogeochemical dynamics for each box:

65. NPD2 � A simple NPD model with two vertical boxes Vertical dynamics:

66. References Blom J.G., Verwer J.G., Journal of Computational and Applied Mathematics, 2000: �A comparison of integration methods for atmospheric transport-chemistry problems� Butensch�n M., 2007: �Numerical Simulations of the coastal marine ecosystem dynamics: integration techniques and data assimilation in a complex physical-biogeochemical model� Haidvogel D. B. and Beckmann A., 1999: �Numerical ocean circulation modelling� Kantha L.H. and Clayson C.A., 2000: �Numerical models of oceans and oceanic processes� Kowalik Z. and Murty T. S., 1993: �Numerical Modelling of Ocean Dynamics� Levy M et al. Geophysical Research Letters, 2001: �Choice of an Advection Scheme for Biogeochemical Models� Madec G., 2008: �NEMO ocean engine� Pietrzak J., Monthly Weather Review, 19998: �The Use of TVD Limiters for Forwart-in-Time Upstream-Biased Advection Schemes in Ocean Modelling� Smolarkiewicz P.K., Numerical Methods in Fluids, 2006: �Multidimensional positive definite advection transport algorithm: An overview�