A Particle tracking schemes
the processes and automatically account for suspension and bed load.
Because of these advantages, Lagrangian schemes have also become more common in the SPM modelling community [Charles et al. , 2008; Krestenitis et al. , 2007; Rolinski et al. , 2005].
Nevertheless most of these models used only small number of particles O(104). Nowadays with easy access to high performance computer clusters, the tracking of individual particles can be parallelised with high efficiency and therefore makes huge particle numbers feasible [Charles et al. , 2008]. This means to deal with particles in the order of >107. This is still negligible, by realising that a bucket of muddy water contains more individual particles, com-pared to the ability of state of the art Lagrangian schemes. Nonetheless increasing the number of particles leads to a better statistical description and makes the answers, a Lagrangian model can give, more reliable.
A.2 The Lagrangian model
holds
hW(t)i= 0 ; Std(W(t)−W(s)) =q|t−s|I (A.4) whereI is the identity matrix. Therefore the noise process has a vanishing mean h·i, its stan-dard deviation scales as √
dt and the increments are uncorrelated.
The first term on the right hand side of Equation (A.3) represents the deterministic part, whereas the second term is the stochastic term. In the case of vanishing turbulent diffusiv-ity, the system of equations reduces to a system of ordinary differential equations (ODEs).
Because the ocean is a turbulent environment, turbulent diffusion has to be included. This is incorporated via the stochastic term. The particles experience a random displacement due to eddies of average size √
2Kdt. Because the turbulent diffusivity K = K(x, t) is spatially highly variable, the term ∇ ·Kneeds to be added to the deterministic part. This corrects for an artificial noise induced drift [Hunter et al. , 1993; Visser, 1997].
A.2.1 Numerical approximation
Because the diffusivity tensor is diagonal, the three spatial directions can be treated separately in developing a numerical approximation to the 3D Langevin equation. Focussing for simplicity on the vertical dimension the following equation needs to be discretised.
dZ(t) = (w + ∂zKZ(z))dt + q2KZ(z) dW(t) (A.5) This equation can further be simplified to
dZ(t) =a(z)dt + b(z)dW(t) (A.6) wherea=w + ∂zKZ(z), represents the deterministic part andb=p2KZ(z) is the stochastic part. Again Eq. (A.6) and (A.5) are only valid in the ˆIto interpretation [Arnold, 1974].
Instead of writing Equation (A.6) in differential form, it is also common to use the integral representation
Zt=Z0 + Z t
0
a(Zs)ds + Z t
0
b(Zs)dWs (A.7)
A straightforward translation of Equation (A.6) into a numerical scheme, is simply to replacedt by ∆tanddW by ∆W. This is equivalent to assuming thata(Zs) andb(Zs) in Equation (A.7) are constant and can be taken out of the integrals. Therefore, the lowest order approximation reads as
Zn+1=Zn + a∆t + b∆Wn (A.8)
This is also known as Euler scheme. In the following, this approximation is named EULER.
This scheme is commonly used [Brickman and Smith, 2001; North et al. , 2006; Spivakovskaya et al. , 2007; Visser, 1997]. Although this is a straightforward approach, some difficulties arise in the
A Particle tracking schemes
case of SDEs. To define the accuracy or order of convergence for stochastic scheme two cases have to be distinguish. For SDEs the order of convergence is separated intoweak and strong [Arnold, 1974; Kloeden and Platen, 1992]. A method is said to have weak/strong order of convergence ofγ if there exists a constant Λ such that
|hp(Zn)i − hp(Z(τ))i| ≤ Λ ∆tγ : weak
h|Zn−Z(τ)|i ≤ Λ ∆tγ : strong (A.9) for any fixed τ =n∆t∈[0, T] and ∆tsufficiently small. Zn represents the true solution and Z(τ) is the approximation. p(·) is an arbitrary function (in most cases a probability density function). The weak criterion asks for the difference in a distribution, whereas the strong criterion accounts for the difference in the trajectory.
As discussed after Equation (A.4) the increment ∆W scales as√
∆t, hence the whole EULER scheme is only of order √
∆t in the strong convergence. Since we are interested in the time evolution of a sediment distribution rather than individual trajectories of single sand grains, the weak convergence is used. In this case the EULER schemes is or order ∆t in the weak sense.
To develop higher order schemes, that have a higher accuracy in the strong definition, the as-sumption thata(Zs) andb(Zs) in Equation (A.7) are constant is not valid any more. Using the appropriate Taylor approximation for the integrals, see e.g. Arnold [1974]; Kloeden and Platen [1992], the next higher order approximation reads as
Zn+1 =Zn + a∆t + b∆Wn + 1
2bb′h(∆Wn)2−∆ti (A.10) where b′ is the spatial derivative. This is also known as the MILSTEIN scheme. This scheme is of order ∆tin the weak and strong convergence. Additional accuracy is gained by including information of the derivative of the noise termb. If this scheme is used in numerical algorithms, the term bb′ can cause problems due to round off errors. To avoid this, a symmetry property of Equation (A.6) can be used. By replacing b again by p2KZ(z) and computing the term bb′ explicitly it turns out that it is equal to∂zKZ(z). After some rewriting, the schemes read as follows
Zn+1 = Zn+w∆t + ∂zK ∆t +√
2K∆W : EULER
Zn+1 = Zn+w∆t + ∂zK (∆W)22+∆t +√
2K∆W : MILSTEIN (A.11)
As stated above the EULER scheme is commonly used, but it is easily appreciated that a higher accuracy is gained here by a simple multiplication with ∆W ·∆W and one addition.
There are no approximations involved and the extra computational cost is negligible.
To further improve the accuracy of the numerical schemes, more terms in the Taylor approxi-mations have to be included [Kloeden and Platen, 1992]. Due to the slow convergence of the
A.2 The Lagrangian model
numerical schemes for SDEs, the extra computational costs are so far prohibitive. Therefore, a multi step scheme is proposed, similar to Runge-Kutta schemes for ODEs.
Zn+1=Zn + 1 2
a( ˜Z) + a∆t + b∆Wn (A.12)
with
Z˜ =Zn + a∆t + b∆Wn (A.13)
Equation (A.12) is a stochastic version of the trapezoidal method also known as HEUN scheme. Note that the predictor step (A.13) is only applied to the deterministic part, the stochastic part is not corrected to keep the numerical approximation consitent with Eq. A.6 [Kloeden and Platen, 1992]. The HEUN scheme, like the MILSTEIN scheme, is of order ∆t in the strong and weak convergence.
At this stage, the question might arise, why three algorithms are presented with the same accuracy. All three algorithms should behave identical in predicting the time evolution of an initial concentration, because they have the same order of convergence. Assume the limit of vanishing diffusivity, here the MILSTEIN scheme becomes identical to the EULER scheme.
This is not the case for the HEUN scheme. Due to the predictor-corrector step, the accuracy is higher. Therefore, in the case of advection-dominated problems, differences will be visible.
This should not be the case if diffusion dominates. However, especially Sec. A.3.3 will show the limitation of the EULER scheme. Moreover, near boundaries, the proper approximations of the particle trajectories become important to resolve for instance the bottom boundary layer.
A.2.2 Boundary conditions
The treatment of boundary conditions is always a critical issue in ocean modelling, especially in coastal regions. The moving sea surface, the sea bottom and lateral boundaries like islands or beaches have to be considered appropriately. E.g., in the framework of PDEs the sea surface is an impermeable boundary and a no flux condition is normally imposed (at least for suspended particulate matter). This no flux condition can be easily violated by overshooting of the trajectories of simulated particles, due to either too large time steps or the random nature of the stochastic increment ∆W. This can lead to a crossing of the boundary. To correct this, a straightforward approach would look like this: When a particle crosses the boundary (due to a too large random displacement), it is simply reflected back into the domain by the amount it penetrates into the boundary domain. It would be advantageous in general to minimise the number of particles that crosses the boundaries in the first place. The first solution that comes to mind is to reduce the time step of the particle displacement. Nevertheless, this would also lead to additional computation time in ’open water’. A more expensive method is the use of a higher-order numerical scheme. This may perhaps not completely prevent the crossing from happening, but it will at least reduce the number of times that it does occur. This was also
A Particle tracking schemes
mentioned by Stijnen et al. [2006].
In the following all boundaries are treated as reflective boundaries (no flux condition), if not stated otherwise.