Numerical Aspects of Marine Biogeochemical Modelling

1
Numerical Aspects of Marine Biogeochemical
Modelling
Wednesday 10th June 2009
Momme Butenschön momm_at_pml.ac.uk
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.

Y(x,t)
Y(p,t)
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.

microscopic
mesoscopic
macroscopic
6
Scales and the validity of the continuum
hypothesis
microscopic
mesoscopic
macroscopic
These scales define the lower limit for the
processes that a continuous model will be able to
resolve.
7
The Balance Equation
Y

• 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

Y
Y
Y
8
The Balance Equation
9
The tracer transport equation
For a tracer in an incompressible fluid, the
terms of efflux and supply in the balance
equation are given by diffusion and local
production/destruction and one obtains
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

forward
backward
central
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.
• gt based on the finite volume formulation of the
  balance equation.

e.g.
fluxes
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-gt0.000000001).

15
Important numerical concepts - Consistency

• For decreasing step-sizes the equations of the
  discrete system approach the original pde.
• Based on local error.
• gt Evaluated on single time step

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 doesnt 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

• Solve

• Single-step schemes
• Euler
• Runge-Kutta

• Multistep schemes
• Adams
• Leap-frog

• Local discretisation error given by
• When expressed as polynomial, lowest order
  polynomial defines order of error

19
The Taylor Theorem

• Given a function and its derivatives are finite
  and continuous, the value at each point tDt is
  given by
• Cut-off error is increasingly small for Dt lt 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
1

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

weight
tn
tn1
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.
• Dont 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.

• 2nd order.
• Steps alternating decoupled, creates
  computational mode that is slowly amplified with
  time.
• gt Needs smoother to converge.

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
The Euler scheme evaluates to
A for the real solution insignificant mode
becomes critical for the discrete solution if the
time step is gt 0.02!
31
Explicit vs. implicit schemes

• Schemes considered so far are explicit, i.e.
  contain only known values on the RHS

• There is in principal no reason, why not to
  include the level n1 on the RHS.
• In the majority of cases this involves more
  complex solution strategies to solve linear or
  even non-linear systems.
• Advantage unconditionally stable.

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.

1
weight
tn
tn1
33
Should I decrease my time step unlimited if I am
not restricted by computational resources?

• No!

discretisation error
rounding error
Error
Dt
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

From Madec G., 2008 "NEMO ocean engine"
41
Central scheme

• 2nd order
• Non-diffusive
• Dispersive
• Unconditionally unstable

From Madec G., 2008 "NEMO ocean engine"
42
Numerical diffusion in upwind scheme

• Applying Taylor series in time and space to the
  upwind scheme leads to

43
Upwind vs. Central Scheme
From Kowalik Z. et al., 1993 Numerical
Modelling of Ocean Dynamics"
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.

From Madec G., 2008 "NEMO ocean engine"
47
Piecewise parallel scheme

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

From Madec G., 2008 "NEMO ocean engine"
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 gt spurious
  transition.
• First order.

52
Operator Splitting Method
Module 1
Module 2
f(c)
f(c)
c
c
I
II
c
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
Module 1
f(c)
I
f
c
c
II
f(c)
Module 2
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?

From GLOBEC, 1993 Sampling and Observational
Systems"
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
• gt 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
• gt 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
P

• Two part exercise
• Spin-up the two boxes to steady state, no
  diffusion.
• Fully coupled dynamics with turbulence cycle.

N
D
P
N
D
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

m
t
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