Barotropic-Baroclinic Mode Splitting

A critical feature of NeuralPOM is the separation of the ocean circulation into fast barotropic (external) modes and slow baroclinic (internal) modes. This mode splitting allows the model to efficiently handle both the rapidly propagating surface gravity waves and the slower three-dimensional density-driven currents.

---

Motivation

The primitive equations support two distinct classes of motion:

Mode	Character	Speed	Time Step Required
Barotropic (external)	Depth-independent, surface gravity waves	$c_0=\sqrt{gH}$ (~200 m/s in deep ocean)	Small ($\mathrm{\Delta }{t}_{\text{BT}}\sim$ 1--10 s)
Baroclinic (internal)	Depth-dependent, density-driven	$\sim$ 1--3 m/s	Larger ($\mathrm{\Delta }{t}_{\text{BC}}\sim$ 100--1000 s)

Without splitting, the entire 3D model would be constrained by the CFL condition of the fastest surface waves. Mode splitting solves the depth-integrated barotropic equations at a much smaller time step than the three-dimensional baroclinic equations, dramatically improving computational efficiency.

---

Depth-Integrated (Barotropic) Equations

The barotropic equations are obtained by vertically integrating the primitive equations from the bottom $z=-H$ to the surface $z=\eta$.

Continuity (Free Surface Evolution)

$$\frac{\partial \eta }{\partial t}+\frac{\partial (D\overline{u})}{\partial x}+\frac{\partial (D\overline{v})}{\partial y}=0$$

where $D=H+\eta$ is the total water depth, and $\overline{u},\overline{v}$ are the depth-averaged velocities:

$$\overline{u}=\frac{1}{D}{\int }_{-H}^{\eta }udz,\overline{v}=\frac{1}{D}{\int }_{-H}^{\eta }vdz$$

The depth-integrated volume transports are:

$$U=D\overline{u}={\int }_{-H}^{\eta }udz,V=D\overline{v}={\int }_{-H}^{\eta }vdz$$

Momentum (Depth-Integrated)

$$\frac{\partial U}{\partial t}+(\text{advection}+\text{diffusion}{)}_U-fV=-gD\frac{\partial \eta }{\partial x}+\frac{{\tau }_s^x-{\tau }_b^x}{{\rho }_0}+F_U^{\text{baroclinic}}$$$$\frac{\partial V}{\partial t}+(\text{advection}+\text{diffusion}{)}_V+fU=-gD\frac{\partial \eta }{\partial y}+\frac{{\tau }_s^y-{\tau }_b^y}{{\rho }_0}+F_V^{\text{baroclinic}}$$

The $F^{\text{baroclinic}}$ terms contain the vertically integrated baroclinic pressure gradient, advection, and horizontal diffusion from the 3D fields, coupling the barotropic and baroclinic modes.

---

Baroclinic (Internal) Equations

The baroclinic equations describe the vertical structure of the flow. They are obtained by subtracting the depth-averaged (barotropic) flow from the full 3D flow:

$$u^{\prime }=u-\overline{u},v^{\prime }=v-\overline{v}$$

The vertical shear evolves according to the 3D momentum equations:

$$\frac{\partial u^{\prime }}{\partial t}+\overrightarrow{v}\cdot \mathrm{\nabla }{u}^{\prime }-f{v}^{\prime }=-\frac{1}{{\rho }_0}\frac{\partial p_{\text{bc}}}{\partial x}+\frac{\partial }{\partial z}\left(K_M\frac{\partial u^{\prime }}{\partial z}\right)+F_u^{\prime }$$

where $p_{\text{bc}}$ is the baroclinic pressure (density-driven component), computed by vertical integration of the horizontal density gradient in baropg().

---

The Free Surface $\eta$ as Coupling Variable

The free surface elevation $\eta (x,y,t)$ is the central variable linking the two modes:

graph TD
    A["3D Baroclinic Fields<br/>(u, v, T, S, ρ)"] -->|"Vertical integration<br/>of baroclinic PGF"| B["Barotropic RHS<br/>(F_U, F_V)"]
    B --> C["Barotropic Solver<br/>(solve for η, U, V)"]
    C -->|"Surface slope<br/>∂η/∂x, ∂η/∂y"| D["3D Momentum<br/>(pressure gradient term)"]
    D -->|"Updated velocities"| A

1. Baroclinic to Barotropic: The vertically integrated baroclinic pressure gradient, advection, and diffusion terms are computed from 3D fields and passed to the barotropic solver as forcing terms.
2. Barotropic to Baroclinic: The updated surface elevation $\eta$ provides the surface slope component of the 3D pressure gradient for the next baroclinic time step.

In the code, the surface pressure gradient in advu is computed as:

# Gradient of Elevation (Future + Back)
grad_eg = (egf[s_i, s_j] - egf[s_im1, s_j]) + (egb[s_i, s_j] - egb[s_im1, s_j])

# Gradient of Atmos Pressure
grad_pa = 2.0 * (e_atmos[s_i, s_j] - e_atmos[s_im1, s_j])

press_grad = grav * 0.125 * _b(dt_sum * (grad_eg + grad_pa) * dy_sum)

Here, egf and egb represent the forward and backward time levels of $\eta$ for the leapfrog integration.

---

Implicit Barotropic Formulation

NeuralPOM's implicit core (implicit_barotropic.py) solves the barotropic equations semi-implicitly, which relaxes the CFL constraint on surface gravity waves.

The Elliptic Equation for $\eta$

The semi-implicit discretization leads to a Helmholtz-type elliptic equation for the free surface:

$$\eta -4\mathrm{\Delta }{t}^{2}g\alpha \mathrm{\nabla }\cdot (D\mathrm{\nabla }\eta )=\text{RHS}$$

In operator form:

$$\mathcal{L}(\eta )=\eta -4\mathrm{\Delta }{t}_i^{2}g\alpha \mathrm{\nabla }\cdot (D\mathrm{\nabla }\eta )$$

where $\alpha$ is the implicitness parameter (typically 1.0 for fully implicit, or 0.5 for Crank-Nicolson).

The elliptic operator is implemented in apply_elliptic_operator:

def apply_elliptic_operator(self, eta, d_2d):
    # H at U/V faces (Arakawa C)
    H_u_east = 0.5 * (d_2d[s_i_p1, s_j] + d_2d[s_i, s_j])
    H_u_west = 0.5 * (d_2d[s_i, s_j] + d_2d[s_i_m1, s_j])
    H_v_north = 0.5 * (d_2d[s_i, s_j_p1] + d_2d[s_i, s_j])
    H_v_south = 0.5 * (d_2d[s_i, s_j] + d_2d[s_i, s_j_m1])

    # Fluxes with land mask (dum @ U-faces, dvm @ V-faces)
    flux_x_east = H_u_east * (eta[s_i_p1, s_j] - eta[s_i, s_j]) / dx_east * dy_east * dum
    flux_x_west = H_u_west * (eta[s_i, s_j] - eta[s_i_m1, s_j]) / dx_west * dy_west * dum
    flux_y_north = H_v_north * (eta[s_i, s_j_p1] - eta[s_i, s_j]) / dy_north * dx_north * dvm
    flux_y_south = H_v_south * (eta[s_i, s_j] - eta[s_i, s_j_m1]) / dy_south * dx_south * dvm

    div_flux = (flux_x_east - flux_x_west + flux_y_north - flux_y_south) / art
    coeff = 4.0 * (self.dti ** 2) * self.grav * self.alpha
    out[s_i, s_j] = eta[s_i, s_j] - coeff * div_flux
    return out

Conjugate Gradient Solver

The elliptic system $\mathcal{L}(\eta )=\text{RHS}$ is solved using a Conjugate Gradient method:

x = (2.0 * self.el - self.elb).clone()  # warm start: linear extrapolation
r = rhs - self.apply_elliptic_operator(x, d_2d)
p = r.clone()
rs_old = (r * r).sum()

for it in range(max_iter):
    Ap = self.apply_elliptic_operator(p, d_2d)
    alpha = rs_old / (p * Ap).sum()
    x += alpha * p
    r -= alpha * Ap
    rs_new = (r * r).sum()
    if rs_new <= tol * tol:
        break
    beta = rs_new / rs_old
    p = r + beta * p
    rs_old = rs_new

A warm start using linear extrapolation in time (${\eta }^{n+1}\approx 2{\eta }^n-{\eta }^{n-1}$) significantly reduces the number of CG iterations required.

Boundary Conditions for $\eta$

The CG solver supports two types of boundary conditions:

Problem Type	BC for $\eta$	Description
iproblem == 4, 5 (Double Gyre)	Dirichlet: $\eta =0$ at boundary	Rigid wall (closed basin)
Others	Neumann: $\frac{\partial \eta }{\partial n}=0$	Ghost cell equals inner cell

def _enforce_eta_neumann_bc(self, eta):
    if self.iproblem == 4 or self.iproblem == 5:
        # Dirichlet: rigid wall
        eta[0, :] = 0.0
        eta[-1, :] = 0.0
        eta[:, 0] = 0.0
        eta[:, -1] = 0.0
    else:
        # Neumann: ghost = inner
        eta[0, :] = eta[1, :]
        eta[-1, :] = eta[-2, :]
        eta[:, 0] = eta[:, 1]
        eta[:, -1] = eta[:, -2]

---

Time Step Ratio

The barotropic time step is typically much smaller than the baroclinic time step:

$$\frac{\mathrm{\Delta }{t}_{\text{barotropic}}}{\mathrm{\Delta }{t}_{\text{baroclinic}}}\equiv \frac{\text{dti}}{\text{dte}}\ll 1$$

Typical values:

Configuration	$\mathrm{\Delta }{t}_{\text{baroclinic}}$	$\mathrm{\Delta }{t}_{\text{barotropic}}$	Ratio
Coastal (1 km, 40 layers)	30--60 s	1--3 s	1:20--1:30
Basin-scale (10 km, 20 layers)	300--600 s	10--20 s	1:30
Deep ocean (25 km, 30 layers)	900--1800 s	30--60 s	1:30

The barotropic solver is called multiple times per baroclinic step (typically 20--30 sub-steps), with the 3D baroclinic forcing terms held constant or linearly extrapolated over the barotropic sub-cycles. In the explicit core, the splitting may use a simpler forward-backward scheme for the barotropic mode.

---

Implementation Architecture

The two cores in NeuralPOM handle the mode splitting differently:

Explicit Core

• Uses a forward-backward time stepping for the barotropic mode
• advave() computes the 2D external mode advection and diffusion at each barotropic sub-step
• The barotropic mode is integrated in time_loop.py with a fixed number of sub-steps per baroclinic step

Implicit Core

• Uses the CG solver for the implicit free surface equation
• solve_implicit_eta_cg() in implicit_barotropic.py handles the elliptic solve
• The $\eta$ solution feeds back into the 3D pressure gradient through egf/egb arrays

The key advantage of the implicit formulation is that the barotropic time step is no longer limited by the CFL condition for surface gravity waves, allowing larger $\mathrm{\Delta }{t}_{\text{barotropic}}$ and fewer sub-steps, or even eliminating sub-stepping entirely.