Boundary conditions

Eastern Boundary: open
The eastern boundary fluxes normal to the faces are specified while values used for advection of momentum and scalars are upwinded

Boundary fluxes
As specified in boundaries.c, the eastern boundary is an open boundary and employs the linearized boundary condition on the velocity with

$$u_b = -h_b\sqrt{\frac{g}{d}} ,.$$

where $h_b$ is the free-surface height, $g$ is gravitational acceleration, and $d$ is the depth, and the minus sign implies flux out of the domain when $h_b>0$. This is enforced in the OpenBoundaryFluxes function with the lines

if(grid->yv[ib]>50) {
   for(k=grid->etop[j];k<grid->Nke[j];k++) 
     ub[j][k]=-phys->h[ib]*sqrt(prop->grav/grid->dv[ib]);
...

The ib index is the index of the cell adjacent to the boundary, and the if-statement is required to distinguish this eastern outflow boundary from the southern boundary. Note that because the edge lengths are $\Delta x = 75$ m, this implies that the distance between the Voronoi points and the boundary edges is $\Delta x/(2\sqrt{3})=21.7$ m. Therefore, we know that if the location of the boundary Voronoi point ib adjacent to a boundary edge j is greater than $50$ m, it must be the eastern boundary.

Boundary velocities
At the eastern open boundary, the boundary velocities which are used for advection of momentum out of the domain are given by the upwind velocity components. This is specified in the BoundaryVelocities function with the lines

if(grid->yv[ib]>50) {
   for(k=grid->etop[j];k<grid->Nke[j];k++) {
      phys->boundary_u[jptr-grid->edgedist[2]][k]=
                                  phys->uc[ib][k];
      phys->boundary_v[jptr-grid->edgedist[2]][k]=
                                  phys->vc[ib][k];
      phys->boundary_w[jptr-grid->edgedist[2]][k]=
                             0.5*(phys->w[ib][k]+phys->w[ib][k+1]);
   }

Because the flux of momentum is calculated at the vertical centers of the vertical faces, the upwind vertical velocity is interpolated from the upper (k) and lower (k+1) faces of the kth cell.

Boundary scalars
The outflow conditions on the scalars are imposed by specified the upwind scalar quantity at the outflow face. This is set in the function BoundaryScalars with the lines

for(k=grid->ctop[ib];k<grid->Nk[ib];k++) {
   phys->boundary_T[jptr-grid->edgedist[2]][k]=phys->T[ib][k];
   phys->boundary_s[jptr-grid->edgedist[2]][k]=phys->s[ib][k];
}

Southern boundary: specified
All quantities at the southern boundary are specified. Since this example simulates a river inflow, the velocity and scalar fields are specified along only part of the southern boundary, while the rest of the boundary is kept at zero inflow to represent a solid boundary. This requires if-statements to determine which of the southern boundary edges fall between the extent of the inflowing river. An alternative is to set the markers of the southern boundary edges that fall within the boundaries of the river to 2, and the rest to 1.

Boundary fluxes Because the velocity components at the southern boundary are specified, the incoming boundary flux is computed in OpenBoundaryFluxes using the velocity components specified in the BoundaryVelocities function, but only within the extent of the 300 m wide river, which is defined for $900 < x < 1200$ m (4 cell faces), as in:

if(grid->xv[ib]>900&&grid->xv[ib]<1200)
  for(k=grid->etop[j];k<grid->Nke[j];k++) 
    ub[j][k]=
      phys->boundary_u[jptr-grid->edgedist[2]][k]*grid->n1[j]+
      phys->boundary_v[jptr-grid->edgedist[2]][k]*grid->n2[j];

Since boundary_u and boundary_v store the inflow velocity components, then this is setting the flux at the inflow to the normal component of the velocity at the inflow, i.e.

$$U_{inflow} = {\mathbf u}_{inflow}\cdot {\mathbf n}_{inflow} ,$$

where ${\mathbf n}_{inflow}$ is, by convention, the inward pointing normal at the boundary face.

Boundary velocities The north-south velocity is specified at the inflow for $900 < x < 1200$ over a specified depth, which represents the river inflow, such that

$$u_{inflow} = 0 ,$$

$$v_{inflow} = \left\{\begin{array}{ll} \mbox{amp} & z>-D_r \\ 0 & \mbox{otherwise ,} \end{array} \right.$$

$$w_{inflow} = 0 .$$

where amp is specified in suntans.dat as $0.01$ m s$^{-1}$, and $D_r=3$ m is the inflow depth. In boundaries.c, this inflow is implemented with

if(grid->xv[ib]>900 && grid->xv[ib]<1200) {
  z=0;
  for(k=grid->etop[j];k<grid->Nke[j];k++) {
    z-=0.5*grid->dzz[ib][k];
    if(z>-3.0)
      phys->boundary_v[jptr-grid->edgedist[2]][k]=prop->amp;
    z-=0.5*grid->dzz[ib][k];
  }
}

All other components are set to zero at the beginning of the main loop in the function BoundaryVelocities with

for(k=grid->etop[j];k<grid->Nke[j];k++) {
  phys->boundary_u[jptr-grid->edgedist[2]][k]=0;
  phys->boundary_v[jptr-grid->edgedist[2]][k]=0;
  phys->boundary_w[jptr-grid->edgedist[2]][k]=0;
}

Boundary scalars The inflow boundary condition for the salinity (or density, when $\beta=1$) field represents an inflow of fresh water with a density anomaly of $\Delta \rho/\rho_0=-10^{-4}$, with a surface depth of $D_r=3$ m, such that

$$\rho_{inflow} = \left\{\begin{array}{ll} \Delta \rho & z>-D_r \\ 0 & \mbox{otherwise .} \end{array} \right.$$

The inflow boundary condition for temperature is a no-gradient condition, and as a result the temperature just outside the boundary is equal to that just inside the boundary. These boundary conditions for salinity (density) and temperature are given in boundaries.c as

if(grid->yv[ib]<50)
  if(grid->xv[ib]>1000 && grid->xv[ib]<1200)
    for(k=grid->ctop[ib];k<grid->Nk[ib];k++) {
      z-=0.5*grid->dzz[ib][k];
      if(z>-3.0) {
        phys->boundary_T[jptr-grid->edgedist[2]][k]=phys->T[ib][k];
        phys->boundary_s[jptr-grid->edgedist[2]][k]=-0.0001;
      } else {
        phys->boundary_T[jptr-grid->edgedist[2]][k]=phys->T[ib][k];
        phys->boundary_s[jptr-grid->edgedist[2]][k]=0;
      }
      z-=0.5*grid->dzz[ib][k];
    }