Numerical implementation of boundary conditions

We adopt a mixed approach for implementing boundary conditions that uses both halo regions and "direct" imposition of boundary conditions, depending on the condition prescribed.

We illustrate how boundary conditions are implemented by considering the tracer equation

where is the diffusive flux of and is an arbitrary source term.

See Model setup: boundary conditions for how to create and use these boundary conditions in Oceananigans.

Gradient boundary conditions

Users impose gradient boundary conditions by prescribing the gradient of a field across an external boundary . The prescribed gradient may be a constant, discrete array of values, or an arbitrary function. The gradient boundary condition is enforced by setting the value of halo points located outside the domain interior such that

where is the vector normal to .

Across the bottom boundary in , for example, this requires that

where   are the heights of the finite volume at and   and  . This prescription implies that the -derivative of across the boundary at   is

as prescribed by the user.

Gradient boundary conditions are represented by the Gradient type.

Value boundary conditions

Users impose value boundary conditions by prescribing , the value of on the external boundary . The value may be a constant, array of discrete values, or an arbitrary function. To enforce a value boundary condition, the gradient associated with the difference between and at boundary-adjacent nodes is diagnosed and used to set the value of the halo point located outside the boundary.

At the bottom boundary in , for example, this means that the gradient of is determined by

which is then used to set the halo point via linear extrapolation.

Value boundary conditions are represented by the Value type.

Flux boundary conditions

Users impose flux boundary conditions by prescribing the flux   of across the external boundary . The flux   may be a constant, array of discrete values, or an arbitrary function. To explain how flux boundary conditions are imposed in Oceananigans.jl, we note that the average of the tracer conservation equation over a finite volume yields

where the surface integral over averages the flux of across the six faces of the finite volume. The right-hand-side of above is denoted as .

An external boundary of a finite volume is associated with a no-penetration condition such that    , where is the vector normal to . Furthermore, the closures currently available in Oceananigans.jl have the property that  . Thus setting     on the external boundary implies that the total flux of across the external boundary is

Oceananigans.jl exploits this fact to define an algorithm that prescribes fluxes across external boundaries :

1. Impose a constant gradient     across external boundaries by using halo points (similar to ), which ensures that the evaluation of in boundary-adjacent cells does not include fluxes across the external boundary, and;

2. Add the prescribed flux to the boundary-adjacent volumes prior to calculating :       , where   denotes values of in boundary-adjacent volumes,   is the flux prescribed along the boundary, is the volume of the boundary-adjacent cell, and is the area of the external boundary of the boundary-adjacent cell. The factor is $-$1 and $+$1 on "left" and "right" boundaries, and accounts for the fact that a positive flux on a left boundary where    implies an "inward" flux of that increases interior values of , whereas a positive flux on a right boundary where   implies an "outward" flux that decreases interior values of .

Flux boundary conditions are represented by the Flux type.

Normal flow boundary conditions

Normal flow boundary conditions directly specify the velocity component normal to the boundary. Typically this is used to impose no penetration boundary conditions, i.e. setting wall normal velocity components to zero on the boundary, or to prescribe the velocity of an inflow or outflow.

On the staggered grid the boundary-normal velocity lies on face points in the boundary direction. For example, the u-component of velocity lies on (Face, Center, Center) points, so on west or east boundaries the point specified by the boundary condition lies on the boundary itself. A tracer on (Center, Center, Center) points has no point on the boundary, and its boundary conditions set the halo points outside of the domain instead.

The other important detail is that normal flow boundary conditions, including no-penetration, are the only conditions used on wall normal velocities when the domain is not periodic. This means that their value affects the pressure calculation for nonhydrostatic models as it is involved in calculating the divergence in the boundary adjacent center point (as described in the fractional step method documentation). Usually boundary points are filled for the predictor velocity (i.e. before the pressure is calculated), and on the corrected field (i.e. after the pressure correction is applied), but for normal flow boundary conditions this would result in the boundary adjacent center point becoming divergent, so they are only filled for the predictor velocity and stay the same after the pressure correction (so the boundary point is filled with the final corrected velocity at the predictor step).

The restriction arises as the boundary condition is specifying the wall normal velocity,  , which leads to the pressure boundary condition

implying that there is a pressure gradient across the boundary. Since we solve the pressure poisson equation (  ) using the method described by , we have to move inhomogeneous boundary conditions on the pressure to the right hand side. In order to do this we define a new field where

This moves the boundary condition to the right hand side as becomes

Given the boundary condition on pressure given above, we can define a new modified predictor velocity which is equal to the predictor velocity within the domain but shares boundary conditions with the corrected field,

The modified pressure Poisson equation becomes   which can easily be solved.

Perhaps a more intuitive way to consider this is to recall that the corrector step projects to the space of divergence free velocity by applying

but we have changed to and to so for   the modified predictor velocity must equal the corrected velocity on the boundary.

For simple normal flow boundary conditions such as no penetration or a straight forward prescription of a known velocity at this is simple to implement as we just set the boundary condition on the predictor velocity and don't change it after the correction. But open boundary conditions, which are normal flow conditions with a scheme, calculate the boundary value based on the interior solution. As a simple example, if we wanted to set the wall normal velocity gradient to zero at the west boundary then we would set the boundary point to

but we then pressure correct the interior so a new error is introduced as

This is preferred to a divergent interior solution as open boundary conditions are typically already unphysical and only used in an attempt to allow information to enter or exit the domain.

Normal flow boundary conditions are represented by the NormalFlow classification.

Open boundary condition "schemes"

Except for trivial cases (i.e. no-penetration) the velocity on the boundary point has to be approximated as it is outside the computed domain. There is insufficient information to step the full equation of motion as gradients across the boundary cannot be computed and simply prescribing a boundary normal velocity is unphysical and reflects energy leaving the domain ().

An open boundary condition is a NormalFlowBoundaryCondition or a ValueBoundaryCondition with a scheme, passed through the scheme keyword argument. The scheme computes the boundary value from the interior solution and is designed to let disturbances propagate out of the domain. In the literature, open boundary conditions of this kind are sometimes called radiation boundary conditions. In a NormalFlowBoundaryCondition the scheme sets the boundary-normal velocity on the boundary face. In a ValueBoundaryCondition it sets the first halo cell of a field that is Center-located in the boundary-normal direction, such as a tracer or a boundary-tangential velocity component. The available schemes are

Scheme	Boundary condition	Fields
PerturbationAdvection	NormalFlowBoundaryCondition, ValueBoundaryCondition	velocities, tracers
NormalRadiation	NormalFlowBoundaryCondition, ValueBoundaryCondition	velocities, tracers
ObliqueRadiation	NormalFlowBoundaryCondition, ValueBoundaryCondition	velocities, tracers
TracerReservoir	ValueBoundaryCondition	tracers
GravityWaveRadiation	NormalFlowBoundaryCondition	barotropic transports U, V
SurfaceWaveRadiation	ValueBoundaryCondition	free surface displacement η

GravityWaveRadiation and SurfaceWaveRadiation apply to the barotropic state of a SplitExplicitFreeSurface.

Perturbation advection

A common method for allowing information to exit is to perform a one-sided advection operation. For example, in the boundary condition fields are advected out of the domain by a locally determined phase speed. We can show that this is the first-order approximation of the full equations of motion in the predictor velocity step. Consider a right boundary normal to the u velocity component (the east boundary):

let    with   where is an externally determined "background" wall-normal flow in the proximity of the boundary, and assume that the divergence of is small,

then, taking only first-order terms:

Now consider the dominant forcing to be relaxation to the background state (we will explain why later) so that   , and recall that when we compute the boundary normal point the interior domain is already stepped, we can take an upwind backward Euler step:

where is the boundary point and and are interpolated to the boundary point. Of course if   (i.e. it is directed into the domain), then we need which is outside of the domain, so this scheme cannot be used. We therefore disregard this term, taking   and  , and rearranging we recover:

For numerical stability, we also need to apply a CFL-like constraint and set  .

Previously we considered the dominant forcing to be relaxation to so that we can take the backward Euler step. And because it is useful to prescribe relaxation on the boundary both to damp numerical oscillations, and to prevent inconsistency when is directed into the domain.

The implementation neglects the terms with and . The relaxation time scale is outflow_timescale when points out of the domain and inflow_timescale when it points into it, and a zero time scale sets the boundary value to . The keyword argument gravity_wave_speed adds a constant outward speed to in , for fields whose signals travel faster than the flow. In a ValueBoundaryCondition the same update applies to a tracer, with computed from the boundary-normal velocity on the boundary face and the tracer relaxed to its exterior value.

Normal radiation

NormalRadiation is the radiation condition of with the adaptive nudging of . A field on the boundary obeys

where is the outward boundary-normal coordinate, is the speed at which signals cross the boundary, is the exterior value given by the boundary condition, and is a relaxation time scale. diagnoses from the solution next to the boundary, where \eqref{eq:normal_radiation} without relaxation gives c_n = - \partial_t \phi / \partial_n \phi.

Let be the value at the boundary point, which is the boundary face in a NormalFlowBoundaryCondition and the first halo cell in a ValueBoundaryCondition, and let and be the values at the first and second interior points next to it. The interior is stepped before the boundary is filled, so and are known when is updated. The phase speed is diagnosed as the Courant number

with   where the denominator vanishes. Inflow and outflow are distinguished by the sign of the boundary-normal velocity . In a NormalFlowBoundaryCondition, is the normal velocity on the first interior face, . In a ValueBoundaryCondition, it is the normal velocity one face into the interior, or on the boundary face when use_boundary_velocity = true. The sign of is not used because is ill-defined when an extremum of crosses the boundary. On outflow the Courant number is bounded below by the advective Courant number and above by one,

and is outflow_timescale. On inflow   and is inflow_timescale. A backward Euler step of on the boundary point gives

with   when  . The update is a convex combination of , and , so it creates no new extrema. With the defaults inflow_timescale = 0 and outflow_timescale = Inf the boundary value is prescribed on inflow and radiates freely on outflow. The adaptive nudging of relaxes strongly on inflow, where the exterior data set the boundary value, and weakly on outflow, where the interior solution does.

The scheme stores and , the values at the start of the time step, at every boundary point. Every stage of a multi-stage time stepper updates from these stored values with the time step of that stage. Before the first time step, is set to , or to where  .

Oblique radiation

ObliqueRadiation extends normal radiation to signals that reach the boundary at an angle, following . With the coordinate tangential to the boundary, the radiation equation becomes

A single equation does not determine both components of the phase velocity. take the phase velocity parallel to , which gives

The Courant numbers are diagnosed at the first interior point from the differences    and   , and from a tangential difference of that is upwinded with respect to the sign of ,

On outflow is bounded as in normal radiation and is bounded to  . On inflow both are zero. The tangential term is discretized explicitly with upwind differences of the boundary values,

where and are backward and forward differences along the boundary. All tangential differences use the values of the previous time step and are zero at the ends of the boundary. Where   the scheme reduces to \eqref{eq:normal_radiation_update}. The tangential term is applied on lateral boundaries only, so top and bottom boundaries use normal radiation. The treatment of inflow and outflow, the advecting velocity, and the stored values are the same as in normal radiation.

Tracer reservoir

A tracer boundary with a prescribed exterior value has no memory of the water that leaves the domain. When the boundary-normal flow reverses, the water that re-enters carries the exterior value regardless of what left through the same face, so oscillating flows such as tides exchange tracer spuriously across the boundary. TracerReservoir follows the tracer reservoirs of MOM6 and keeps a reservoir value just outside each boundary point. The reservoir relaxes toward an upstream value over a distance advected across the boundary,

where is the boundary-normal velocity on the boundary face. On outflow  , the adjacent interior value, and is outflow_length_scale. On inflow   and is inflow_length_scale. A backward Euler step gives

which is stable for any , and the halo cell is set to . Water that re-enters the domain therefore carries back the tracer that recently left, and relaxes to over an advected distance of about inflow_length_scale. A zero length scale sets   at every time step and an infinite one freezes the reservoir. With both length scales zero, the default, TracerReservoir reduces to a condition without memory: zero gradient on outflow and on inflow. The reservoir starts from and is advanced over the stages of a time step in the same way as the stored values of normal radiation.

Gravity wave radiation

A SplitExplicitFreeSurface advances the barotropic transports    and    and the free surface displacement with many short substeps of the barotropic equations. Their open boundaries must let surface gravity waves, which travel at   in a water column of depth , leave the domain, and let exterior signals such as tides enter it. Linearized in the direction normal to an east boundary, the barotropic equations are

whose Riemann invariants   and   travel at the speeds and . At an east boundary the first leaves the domain and the second enters it. GravityWaveRadiation is the condition of , which takes the incoming invariant from the exterior state,    , so that

On west and south boundaries the incoming invariant is   and the sign of the second term is reversed. The exterior transport and free surface are the 2-tuple (U, η) value of the boundary condition, see GravityWaveRadiationBoundaryCondition. The free surface on the boundary face, , is the average of the first interior and halo values of , is the column depth on the boundary face, and   where the column is dry.

The condition is applied at every barotropic substep. A target_transport shifts the transports on the boundary face after every application of the condition,

where is the length element along the boundary and is the length of its wet part. After the substeps, the barotropic correction adjusts the three-dimensional velocity on the boundary face so that its depth integral equals the barotropic transport. The boundary condition on the three-dimensional normal velocity, for example a NormalRadiation scheme, therefore sets only the departure from the depth average.

Surface wave radiation

SurfaceWaveRadiation is the companion condition on , from . The free surface radiates at the gravity wave speed,

which is discretized implicitly as in ROMS,

where is the halo value, the adjacent interior value, the column depth of the interior cell and the barotropic substep. The halo value evolves with the radiation equation instead of mirroring the interior value, as the default zero-gradient condition on does, and it enters the condition through the face average .

In each barotropic substep the halos are filled first, then and are advanced and their boundary values set by GravityWaveRadiation, and finally is advanced with the new transports. A SplitExplicitFreeSurface adds SurfaceWaveRadiation on at every boundary where or has a GravityWaveRadiation condition, unless boundary conditions for are given. tidal_boundary_conditions builds both conditions on every lateral boundary from the harmonic constants of a tidal atlas.

Net transport through open boundaries

Before the pressure solve of a NonhydrostaticModel, the boundary-normal velocities set by the schemes are corrected so that the net volume transport into the domain vanishes. The target_transport keyword argument of the schemes gives the user more control on how the correction is distributed among the boundaries.

The correction has two steps. First, the normal velocity on every boundary whose scheme has a target_transport is shifted uniformly so that the transport through that boundary matches the target,

where is the area of the boundary and the transport is positive in the positive coordinate direction. Second, the remaining net inflow is removed by the same uniform shift of the normal velocity on all boundaries that have a scheme and no target. If every boundary with a scheme has a target, the targets and the transport through boundaries with a prescribed normal velocity must add up to zero net transport into the domain.

A HydrostaticFreeSurfaceModel skips the second step because the free surface absorbs a net imbalance. A net transport through the open boundaries changes the mean free surface height, so the targets do not need to add up to zero net transport. With a SplitExplicitFreeSurface the transport can be controlled through the target_transport of GravityWaveRadiation on the barotropic transport.