Time-\/averaged accelerations

The barotropic equations are timestepped with a timestep to resolve the surface gravity waves. With care, the baroclinic timestep only need resolve the inertial oscillations. The barotropic accelerations are averaged over the many barotropic timesteps taken between baroclinic steps. At time $n$, the baroclinic accelerations are computed. The vertical average of that acceleration is subtracted off and replaced by the time-\/averaged acceleration from the group of barotropic timesteps\+:

\[ \Delta t \frac{\partial \vec{u}}{\partial t} = \Delta t \left( \frac{\partial \vec{u}}{\partial t} - \frac{\partial \vec{u}_{BT}}{\partial t} \right)^n + \Delta t \overline{\frac{\partial \vec{u}_{BT}}{\partial t}}^{\Delta t} \]

Similarly, the velocities used in the tracer equation are a careful blend of the barotropic and baroclinic solutions\+:

\[ \Delta t \frac{\partial \theta}{\partial t} + \Delta t \left( \tilde{\vec{u}} \cdot \nabla \theta + \widetilde{w} \frac{\partial \theta}{\partial z} \right) \] with \[ \tilde{\vec{u}} = \vec{u}_{BC} + \overline{\vec{u}_{BT}}^{\Delta t} \] \[ \frac{\partial \widetilde{w}}{\partial z} = - \nabla \cdot \tilde{\vec{u}} \]\hypertarget{Barotropic_Baroclinic_Coupling_SSH_Estimates}{}\doxysection{Two estimates of the free surface height}\label{Barotropic_Baroclinic_Coupling_SSH_Estimates}
The vertically discrete, temporally continuous layer continuity equations are\+:

\[ \frac{\partial h_k}{\partial t} = - \nabla \cdot (\vec{u} h_k) = - \nabla \cdot \mathbf{F} (u_k, h_k) \]

The relationship between the free surface height and the layer thicknesses $h_k$ is\+:

\[ \eta = \sum_{k=1}^N h_k - D \]

We get an evolution equation for the free surface height\+:

\[ \frac{\partial \eta}{\partial t} = \sum_{k=1}^N \frac{\partial h_k}{\partial t} = - \nabla \cdot \sum_{k=1}^N \mathbf{F} (u_k, h_k) \]

If the algorithms for the fluxes in the continuity equations are {\itshape linear} in the velocity, the free surface height can be rewritten as\+:

\[ \frac{\partial \eta}{\partial t} \&= - \nabla \cdot \sum_{k=1}^N \mathbf{F} (u_k, h_k) = - \nabla \cdot \sum_{k=1}^N (\vec{u}_k h_k) \\ \&= - \nabla \cdot \left[ \sum_{k=1}^N h_k \frac{\sum_{k=1}^N (\vec{u}_k h_k)} {\sum_{k=1}^M k_k} \right] \equiv - \nabla \cdot H \mathbf{U} \] where

\[ \mathbf{U} \equiv \frac{\sum_{k=1}^N (\vec{u}_k h_k)} {\sum_{k=1}^M k_k} \] \[ H \equiv \sum_{k=1}^N h_k \] However, A\+LE models like M\+O\+M6 require positive-\/definite, nonlinear continuity solvers (M\+O\+M6 uses \mbox{\hyperlink{PPM}{P\+PM Advection Scheme}}). We need a different way to reconcile this estimate of free surface height with the one coming from the barotropic equations. After rejecting several other options, M\+O\+M6 is now using the scheme from \cite{hallberg2009}. The barotropic update of $\eta$ is given by\+:

\[ \frac{\eta^{n+1} - \eta^n}{\Delta t} + \nabla \cdot \left( \overline{UH} \right) = 0 \]

Determine the $\Delta U$ such that\+:

\[ \sum_{k=1}^N \mathbf{F} (\tilde{u}_k, h_k) = \left( \overline{UH} \right) \] where \[ \tilde{u}_k = u_k + \Delta U \]

The layer timestep equation becomes\+:

\[ h_k^{n+1} = h_k^n - \Delta t \nabla \cdot \mathbf{F} (\tilde{u}_k, h_k) \]

This scheme has these properties\+:

\begin{DoxyItemize}
\item Total mass is conserved layer-\/wise. \item Layer equations retain their flux form. \item Total salt, heat, and other tracers are exactly conserved. \item Free surface heights exactly agree. \item Requires (very few) completely local iterations. \item The velocity corrections are barotropic, and more likely to correspond to the layers whose flow was deficient than in some older schemes.\end{DoxyItemize}
The solution is unique provided that\+: \[ \frac{\partial}{\partial \tilde{u}_k} \mathbf{F}(\tilde{u}_k, h_k) > 0 \] This is a reasonable requirement for any discretization of the continuity equation. In the continuous limit, $\mathbf{F} (u,h) = uh$, so one interpretation is\+: \[ \frac{\partial}{\partial \tilde{u}_k} \mathbf{F}(\tilde{u}_k, h_k) = h_{k,\mbox{Marginal}} \] With the P\+PM continuity scheme\+: \[ F_{i+1/2} = \frac{1}{\Delta t} \int_{x_{i+1/2} - u \Delta t}^{x_{i+1/2}} h_i^n (x) dx \] leads to\+: \[ \frac{\partial F_{i+1/2}}{\partial u_{i+1/2}} = h_i^n ( x_{i+1/2} - u_{i+1/2} \Delta t ) \equiv h_{k, \mbox{Marginal}} \] $h_i(x) > 0$ is already required for positive definiteness.

Newton\textquotesingle{}s method iterations quickly give $\Delta U$\+: \[ \Delta U^{m+1} = \Delta U^m + \frac{\left( \overline{UH} \right) - \sum_{k=1}^N F(u_k + \Delta U^m, h_k)}{\sum_{k=1}^N h_{k,\mbox{Marginal}}} \]\hypertarget{Barotropic_Baroclinic_Coupling_subsec_practical}{}\doxysubsection{How practical is this iterative approach?}\label{Barotropic_Baroclinic_Coupling_subsec_practical}
The piecewise parabolic method continuity solver uses two steps\+:

\begin{DoxyItemize}
\item Set up the positive-\/definite subgridscale profiles, $h_{PPM}(x)$.\end{DoxyItemize}
 \begin{DoxyImage}\includegraphics[width=\textwidth,height=\textheight/2,keepaspectratio=true]{h_PPM.png}\doxyfigcaption{Piecewise parabolic reconstructions of $h(x)$.}\end{DoxyImage}

\begin{DoxyItemize}
\item Integrating the profiles to determine $F$. Step 1 is typically $\sim 3$ times as expensive as step 2. $F(u,h)$ is piecewise cubic in $u$, but often nearly linear. Newton\textquotesingle{}s method iterations converge quickly\+:\end{DoxyItemize}
 \begin{DoxyImage}\includegraphics[width=\textwidth,height=\textheight/2,keepaspectratio=true]{Newton_PPM.png}\doxyfigcaption{Newton's method iterations for finding $\Delta U$.}\end{DoxyImage}

In a global model the sea surface heights converge everywhere to a tolerance of $10^{-6}$ m within five iterations. These five iterations add $\sim 1.6$ times more C\+PU time to the P\+PM continuity solver and the continuity solver is just 12\% of the total model time.\hypertarget{Barotropic_Baroclinic_Coupling_bottom_drag}{}\doxysubsection{A note on bottom drag}\label{Barotropic_Baroclinic_Coupling_bottom_drag}
Barotropic accelerations do not lead to barotropic flows after a timestep when vertical viscosity is taken into account.

\[ u_k^{n+1} = u_k^n + \Delta t A_k + \Delta t \frac{\tau_{k-1/2} - \tau_{k+1/2}}{h_k} \]

\[ \tau_{1/2} = \tau_{Wind} \] \[ \tau_{k+1/2} = \nu_{k+1/2} \frac{u_k^{n+1} - u_{k+1}^{n+1}}{h_{k+1/2}} \] \[ \tau_{N+1/2} = \nu_{N+1/2} \frac{2 u_N^{n+1}}{h_{N+1/2}} \] \[ \gamma_k \equiv \frac{1}{\Delta t} \frac{\partial u_k^{n+1}} {\overline{\partial A}} \]

A tridiagonal equation for $\gamma_k$ results, going from 0 for thin layers near the bottom to 1 far above the bottom.

\[ \gamma_1 = 1 + \frac{1}{h_1} \left[ - \frac{\nu_{3/2} \Delta t}{h_{3/2}} (\gamma_1 - \gamma_2) \right] \] \[ \gamma_k = 1 + \frac{1}{h_k} \left[ \frac{\nu_{k-1/2} \Delta t}{h_{k-1/2}} (\gamma_{k-1} - \gamma_k) - \frac{\nu_{k+1/2} \Delta t}{h_{k+1/2}} (\gamma_k - \gamma_{k+1}) \right] \] \[ \gamma_N = 1 + \frac{1}{h_N} \left[ \frac{\nu_{N-1/2} \Delta t}{h_{N-1/2}} (\gamma_{N-1} - \gamma_N) - \frac{2\nu_{N+1/2} \Delta t}{h_{N+1/2}} \gamma_N \right] \]

In the continuous limit\+:

\[ \gamma(z) = 1 + \Delta t \frac{d}{dz} \left( \nu \frac{d \gamma}{dz} \right) \] with boundary conditions\+: \[ \frac{d \gamma}{dz} (0) = 0 \] \[ \gamma(-D) = 0 \]

For deep water and constant $\nu$, we get\+: \[ \gamma (z) = 1 - \exp \left( - \sqrt{\nu \Delta t} (z + D) \right) \]

 
\begin{DoxyImage}
\includegraphics[width=\textwidth,height=\textheight/2,keepaspectratio=true]{bottom_drag.png}
\doxyfigcaption{The continuous solution for barotropic flow plus a no-\/slip condition at the bottom.}
\end{DoxyImage}


We want a scheme in which the split model looks exactly like the unsplit model would if we had taken all those short 3D timesteps. Rather than applying a barotropic velocity, we use a barotropic acceleration and allow the continuity solver to determine the transports consistent with a no-\/slip bottom boundary layer and perhaps also a no-\/slip surface boundary layer under an ice shelf.

From above, the barotropic equation is\+:

\[ \frac{\eta^{n+1} - \eta^n}{\Delta t} + \nabla \cdot \left( \overline{UH} \right) = 0 \] We need to determine the $\Delta \overline{A}$ such that\+: \[ \sum_{k=1}^N \mathbf{F} (\tilde{u}_k, h_k) = \left( \overline{UH} \right) \] where \[ \tilde{u}_k = u_k + \gamma_k \Delta \overline{A} \Delta t \]\hypertarget{Barotropic_Baroclinic_Coupling_bt-bc_details}{}\doxysection{Additional details about the split time stepping}\label{Barotropic_Baroclinic_Coupling_bt-bc_details}
\begin{DoxyItemize}
\item Transports are used as input and output to the barotropic solver. The continuity solver is inverted to determine velocities.\end{DoxyItemize}
\[ \frac{\partial \eta}{\partial t} = \nabla \cdot \overline{U} + M \] \[ \overline{U} (\overline{u}) = \frac{1}{\Delta t} \int_0^{\overline{u} \Delta t} H(x) dx \] \[ \overline{u}^n = \overline{U}^{-1} \left( \sum_{k=1}^N U_k^n \right) \] \[ u_k^{n+1} = \tilde{u}_k^{n+1} + \Delta \overline{u} \]

We need to find $\Delta \overline{u}$ such that\+:

\[ \sum_{k=1}^N U_k \left( \tilde{u}_k^{n+1} + \Delta \overline{u} \right) = \overline{U}^{n+1} \]

\begin{DoxyItemize}
\item Barotropic accelerations are treated as anomalies from the baroclinic state\+:\end{DoxyItemize}
\[ \frac{\partial \overline{u}}{\partial t} \&= - f \hat{k} \times (\overline{u} - \overline{u}_{Cor}) - \nabla \overline{g} (\eta - \eta_{PF}) - \\ \& \frac{c_D \left( ||u_{Bot}|| + ||u_{Shelf}|| \right)}{\sum_{k=1}^N h_k} (\overline{u} - \overline{u}_{Drag}) + \frac{\sum_{k=1}^N h_k \frac{\partial u_k}{\partial t}} {\sum_{k=1}^N h_k} \]

\begin{DoxyItemize}
\item Bottom drag (and under ice surface-\/drag) is treated implicitly. \item The barotropic continuity solver uses flow-\/dependent thickness fits which approximate the sum of the layer thickness transports, to accommodate wetting and drying. An image of this is shown here\+:\end{DoxyItemize}
 
\begin{DoxyImage}
\includegraphics[width=\textwidth,height=\textheight/2,keepaspectratio=true]{bt_bc_thickness.png}
\doxyfigcaption{The barotropic transports depend on the baroclinic flows and thicknesses.}
\end{DoxyImage}
\hypertarget{Barotropic_Baroclinic_Coupling_time-split_summary}{}\doxysection{Summary of M\+O\+M6 split time stepping}\label{Barotropic_Baroclinic_Coupling_time-split_summary}
\begin{DoxyItemize}
\item We use an efficient approach for handling fast modes via simplified 2-\/d equations, while the 3-\/d baroclinic timestep is determined by baroclinic dynamics.\end{DoxyItemize}
\begin{DoxyItemize}
\item The barotropic solver determines free surface height and time-\/averaged depth-\/integrated transports.\end{DoxyItemize}
\begin{DoxyItemize}
\item By using anomalies, the M\+O\+M6 split solver gives identical answers to an equivalent unsplit scheme for short timesteps.\end{DoxyItemize}
\begin{DoxyItemize}
\item This scheme has been demonstrated to work with wetting and drying, as well as under ice-\/shelves.\end{DoxyItemize}
\begin{DoxyItemize}
\item The linear barotropic solver allows M\+O\+M6 to automatically set a stable barotropic timestep (e.\+g. to 98\% of maximum). \end{DoxyItemize}