Numerical methods in OceanTurb.jl
OceanTurb.jl
uses a one-dimensional finite-volume method to discretize momentum, temperature, salinity, and other variables in the $z$-direction. A variety of explicit and implicit-explicit schemes are implemented for temporal integration.
Spatial discretization
An ASCII-art respresentation of an example grid with N=3
is
▲ z
|
i=4 *
j=4 === Top ▲
i=3 * | Δf[3]
j=3 --- ▼
i=2 * ▲
j=2 --- | Δc[2]
i=1 * ▼
j=1 === Bottom
i=0 *
where the double lines indicate the top and bottom of the domain, the single lines indicate "face" boundaries, the i
's index cell centers (nodes) and j
's index the $z$-location of cell interfaces (faces). Horizontal momentum and tracer variables are located at cell centers, while fluxes of these quantities (and vertical-velocity-like variables when present) are located at cell faces. The cells at $i=0$ and $i=4$ are 'ghost cells', whose values are set according to the boundary condition. For a no flux or zero gradient boundary condition, for example, we would set c[0]=c[1]
and c[4]=c[3]
.
Finite volume derivatives and fluxes
The derivative of a quantity $\Phi$ at face $i$ is
where $\Phi_i$ denotes the value of $\Phi$ at cell $i$, and $\Delta c_i = z_{c, i} - z_{c, i-1}$ is the vertical separation between node $i$ and cell point $i-1$.
With diffusivities defined on cell interfaces, the diffusive flux across face $i$ is $K_i \left ( \d_z \Phi \right )_i$. The (negative of the) divergence of the diffusive flux at node $i$ is therefore
In the top cell where $i=N$, the diffusive flux is
In the bottom cell where $i=1$, on the other hand, the diffusive flux is
With negative advective velocities $A$ (corresponding to downdrafts or down-travelling plumes), defined at cell centers, the advective flux divergence is
Time integration
To integrate ocean surface boundary layer models forward in time, we implement various explicit and implicit-explicit time-stepping schemes. The function time_step!(model, Δt, Nt)
steps a model forward in time.
Timesteppers in OceanTurb.jl
integrate equations of the form
where $\Phi(z, t)$ is a variable like velocity, temperature, or salinity, $A$, $K$, and $L$ are an advective 'mass flux', diffusivity, and damping coefficient which are a general nonlinear functions of $\Phi$, $z$, and external parameters, and $R$ is an arbitrary function representing any number of processes, including the Coriolis force or external forcing.
Time integration methods
We implement time_step!
functions and types for:
- explicit forward Euler
- semi-implicit backward Euler
Forward Euler method
The explicit forward Euler time integration scheme obtains $\Phi$ at time-step $n+1$ using the formula
where the superscripts $n$ and $n+1$ denote the solution at time-step $n$ and $n+1$, respectively.
Backward Euler method
The backward Euler method obtains $\Phi$ at time-step $n+1$ using the formula
The $z$-derivatives in the advective and diffusive terms generate an elliptic problem to be solved for $\Phi^{n+1}$ at each time-step. In the finite volume discretization used by OceanTurb.jl
, this elliptic problem becomes a matrix problem of the form
where $\mathcal{L}^n_{ij}$ is a matrix operator at time-step $n$, and the subscripts $i$ or $j$ denote grid points $i$ or $j$. For the diffusive problems considered by our backward Euler solver, the matrix multiplication $\mathcal{L}^n_{ij} \Phi_j^{n+1}$ has the form
To form the matrix operator in \eqref{implicitoperatormatrix}, we have used the second-order flux divergence finite difference operators in \eqref{fluxdivop}–\eqref{fluxdivop_bottom}.
It is crucial to note that the diffusive operator that contributes to $\mathcal{L}^n_{ij}$ does not include fluxes across boundary faces. In particular, $\mathcal{L}^n_{ij}$ in \eqref{implicitoperatormatrix} enforces a no-flux condition across the top and bottom faces. Accordingly, fluxes through boundary faces due either to Dirichlet (Value) boundary conditions or non-zero fluxes are accounted for by adding the contribution of the flux diverence across the top and bottom face to $R \left ( \Phi \right )$.