Shallow Water Model Formulation

J. Emmanuel Johnson
CNRS
MEOM

Shallow Water

Shallow Water Equations PDE.

Figure 1:The shallow water equations schematic as per [Vallis, 2019].

h(x,y,t)=η(x,y,t)ηB(x,y,t)h(x,y,t) = \eta(x,y,t) - \eta_B(x,y,t)

where hh is the total fluid thickness, η\eta is the height of the upper free surface, and ηB\eta_B is the height of the lower surface or the bottom topography. Note: when the lower surface is completely flat, then ηB\eta_B will be equivalent to zero so h=ηh=\eta exactly. However, when we do have bottom topography or a free form layer beneath, η2\eta_2, we need to include the ηB\eta_B

We can also define the layer thickness to be a function of the surface displacement and mean height as

h(x,y,t)=η(x,y,t)+H(x,y)h(x,y,t) = \eta(x,y,t) + H(x,y)

where HH is the mean thickness. Note, all methods are related to one another via this formula

η=h+ηB=H+ηT\eta = h + \eta_B = H + \eta_T

where ηT\eta_T is the deviation from the top of the total mean thickness.


NonLinear Form

For the SW equations, we define a PDE that captures spatiotemporal relationships between the total thickness height, hh, and well as the velocity components (u,v)(u,v).

Height:h=h(x,t)xΩRDs,tTR+Zonal Velocity:u=u(x,t)xΩRDs,tTR+Meridonal Velocity:v=v(x,t)xΩRDs,tTR+\begin{aligned} \text{Height}: && h & =\boldsymbol{h}(\vec{\mathbf{x}},t) && && \vec{\mathbf{x}}\in\Omega\sub\mathbb{R}^{D_s} , && t\in\mathcal{T}\sub\mathbb{R}^{+} \\ \text{Zonal Velocity}: && u & =\boldsymbol{u}(\vec{\mathbf{x}},t) && && \vec{\mathbf{x}}\in\Omega\sub\mathbb{R}^{D_s} , && t\in\mathcal{T}\sub\mathbb{R}^{+} \\ \text{Meridonal Velocity}: && v & =\boldsymbol{v}(\vec{\mathbf{x}},t) && && \vec{\mathbf{x}}\in\Omega\sub\mathbb{R}^{D_s}, && t\in\mathcal{T}\sub\mathbb{R}^{+} \\ \end{aligned}

We can write out the momentum constraints for the velocity vector-field, u\vec{\boldsymbol{u}}, and the height field, hh, explicitly. For a single-layer fluid, including the Coriolis term, and the inviscid SW equations as

Momentum:Dtu+f×u=g(h+ηB)Mass:th+(hu)=0\begin{aligned} \text{Momentum}: && && D_t\vec{\boldsymbol{u}} + f\times \vec{\boldsymbol{u}} &= - g \boldsymbol{\nabla}(h + \eta_B) \\ \text{Mass}: && && \partial_t h + \boldsymbol{\nabla}\cdot(h\vec{\boldsymbol{u}}) &= 0 \end{aligned}

where huh\vec{\boldsymbol{u}} is the horizontal velocity, ff is the planetary vorticty (equation (28)) and gg is the gravitational constant. We can expand these equations to be:

Height:th+x(hu)+y(hv)=0Zonal Velocity:tu+uxu+vyufv=gx(h+ηB)+Fx+Bx+Mx+ξxMeridonal Velocity:tv+uxv+vyv+fu=gy(h+ηB)+Fy+By+My+ξy\begin{aligned} \text{Height}: && && \partial_t h + \partial_x (hu) + \partial_y (hv) &= 0\\ \text{Zonal Velocity}: && && \partial_t u + u \partial_x u + v\partial_y u - fv &= - g \partial_x (h + \eta_B) + F_x + B_x + M_x + \xi_x \\ \text{Meridonal Velocity}: && && \partial_t v + u \partial_x v + v\partial_y v + fu &= - g \partial_y (h + \eta_B) + F_y + B_y + M_y + \xi_y \end{aligned}

We have a lot of extra terms in the above equation as was proposed in [Klöwer et al., 2018]. These take into account all of the external and/or internal forces and dissipations. We outline each of them explicitly below.

Forcing Terms

There is a forcing vector which operates on the (x,y)(x,y) plane and changes wrt time, i.e. f:Ω×TR2\vec{\boldsymbol{f}}:\Omega\times\mathcal{T} \rightarrow \mathbb{R}^{2}. Each component of the vector operates on an individual velocity component as shown in (6). In general, it is defined as

f=[Fx,Fy]=[Fx(x,t),Fy(x,t)]\vec{\boldsymbol{f}} = [F_x, F_y]^\top = \left[ \boldsymbol{F}_x(\vec{\boldsymbol{x}},t), \boldsymbol{F}_y(\vec{\boldsymbol{x}},t) \right]^\top

This forcing operates on the top of the model, for example the wind or atmospheric component.

Bottom Friction Term

b=[Bx,By]=[Bx(u,h),By(u,h)]\vec{\boldsymbol{b}} = [B_x, B_y]^\top = \left[ \boldsymbol{B}_x(\vec{\boldsymbol{u}},h), \boldsymbol{B}_y(\vec{\boldsymbol{u}},h) \right]^\top

Lateral Mixing of Momentum

m=[Mx,My]=[Mx(u,h),My(u,h)]\vec{\boldsymbol{m}} = [M_x, M_y]^\top = \left[ \boldsymbol{M}_x(\vec{\boldsymbol{u}},h), \boldsymbol{M}_y(\vec{\boldsymbol{u}},h) \right]^\top

Negative Viscosity Backscatter

Ξ=[ξx,ξy]\vec{\boldsymbol{\Xi}} = [\xi_x, \xi_y]^\top

Vector Invariant Formulation

We can rewrite the formulation for the momentum parts, (u,v)(u,v), in equation (6) in the vectorized notation

tu+uu+f×u=gh\partial_t \vec{\boldsymbol{u}} + \vec{\boldsymbol{u}} \cdot \boldsymbol{\nabla} \vec{\boldsymbol{u}} + f\times \vec{\boldsymbol{u}} = - g \boldsymbol{\nabla}h

We will use the vector identity of equation (33) and plug this into equation (11) to arrive at the vector form.

Height:th=x(hu)y(hv)=0Zonal Velocity:tu=qhvxp+Fx+Mx+Bx+ξxMeridonal Velocity:tv=qhuyp+Fy+My+By+ξy\begin{aligned} \text{Height}: && && \partial_t h &= -\partial_x (hu) - \partial_y (hv) = 0\\ \text{Zonal Velocity}: && && \partial_t u &= qhv - \partial_x p + F_x + M_x + B_x + \xi_x \\ \text{Meridonal Velocity}: && && \partial_t v &= - qhu - \partial_y p + F_y + M_y + B_y + \xi_y \end{aligned}

where qq is the potential vorticity (equation (30)) and pp is the work given by the Bernoulli potential (equation (32)).


Linearized Shallow Water Equations

We can remove the advection terms from equation (6) by linearizing the PDE.

th+H(xu+yv)=0tufv=gxhκutv+fu=gyhκv\begin{aligned} \partial_t h &+ H \left(\partial_x u + \partial_y v \right) = 0 \\ \partial_t u &- fv = - g \partial_x h - \kappa u \\ \partial_t v &+ fu = - g \partial_y h - \kappa v \end{aligned}

This is appropriate for when the Rossby number is small. We also assume that the wave height is much smaller than the actual mean height.

Example Applications

Data Assimilation

This has been used in data assimilation schemes [Guillou et al., 2021] to jointly assimilate observations of SSH whereby the a QG model was used for the balanced motions and a linear SW mode was used for the internal tides.


Wind Forcing

We can construct a very typical wind forcing system that is typical within the Northern Hemisphere at mid-latitudes. It will resemble the tradewinds in the southern part of the domain and the westerlies in the northern part.

We first define a domain with a height, HH and a basin length, (Lx,Ly)(L_x,L_y).

Fy=0Fx=F0ρ0H[cos(2π(2Ly12))+2sin(2π(yLy12))]\begin{aligned} F_y &= 0 \\ F_x &= \frac{F_0}{\rho_0 H} \left[ \cos \left( 2\pi \left(\frac{2}{L_y} - \frac{1}{2} \right)\right) + 2\sin \left(2\pi \left( \frac{y}{L_y} - \frac{1}{2} \right) \right) \right] \end{aligned}

where F0=0.12F_0=0.12Pa and ρ0=1e3\rho_0=1e3 kgm3^{-3}.


Boundary Conditions

Kinematic Boundaries. These boundaries prevent and flow from going through the boundaries. These are typically on the East-West extent as the flow moves from left/right to right/left

u(x=0)=u(x=Lx)=0v(x=0)=v(x=Ly)=0\begin{aligned} u(x=0)=u(x=L_x)&=0 \\ v(x=0)=v(x=L_y)&=0 \end{aligned}

Gradient Boundaries. We can also say that there is no gradient across the boundaries for the height, hh. So we can define this as

xh(x=0)=xh(x=Lx)=0yh(x=0)=yh(x=Ly)=0\begin{aligned} \partial_x h(x=0)=\partial_x h(x=L_x)&=0 \\ \partial_y h(x=0)=\partial_y h(x=L_y)&=0 \end{aligned}

No-Slip Boundaries. These refer to the idea that the tangential velocity should vanish. We typically prescribe these for the North-South extent.

u(y=0)=u(y=Ly)=0v(y=0)=v(y=Lx)=0\begin{aligned} u(y=0)=u(y=L_y)&=0 \\ v(y=0)=v(y=L_x)&=0 \end{aligned}

Free-Slip Boundaries. These refer to the idea that the tangential velocity gradients should be zero. We often refer to these as friction-less boundaries. We also typically prescribe these on the North-South extent.

yu(y=0)=yu(y=Ly)=0xv(x=0)=xv(x=Lx)=0\begin{aligned} \partial_y u(y=0)=\partial_y u(y=L_y)&=0 \\ \partial_x v(x=0)=\partial_x v(x=L_x)&=0 \end{aligned}

HelpFul Tools

Material Derivative

This is an operator that models the movement of a fluid parcel within a Eulerian framework.

Dt=t+u=t+ux+vyD_t = \partial_t + \vec{\boldsymbol{u}} \cdot \boldsymbol{\nabla} = \partial_t + u\partial_x + v\partial_y

Planetary Vorticity

f(y)=2Ωsin(θ)+1R2Ωcos(θ)yf(y) = 2\Omega\sin(\theta) + \frac{1}{R}2\Omega\cos(\theta) y

where f0=2Ωsin(θ0)f_0=2\Omega\sin(\theta_0) is the Coriolis force and β=1R2Ωcos(θ0)\beta=\frac{1}{R}2\Omega\cos(\theta_0) is the approximate β\beta-plane forcing. For more information, see this wiki article for a better overview or this video for a more in-depth overview.

Relative Vorticity

ζ=xvyu=×u\zeta = \partial_x v - \partial_y u = \boldsymbol{\nabla}\times\vec{\boldsymbol{u}}

Potential Vorticity

q=ζ+fhq = \frac{\zeta + f}{h}

where ζ\zeta is the relative vorticity given by (29), ff is the planetary vorticity given by (28). A good overview between the relationshop between each of the vorticity equations can be found in this youtube video.

Kinetic Energy

KE=12(u2+v2)\mathcal{KE} = \frac{1}{2}\left( u^2 + v^2\right)

Bernoulli Potential

Bernoulli Potential:p=12(u2+v2)+g(h+ηB)=KE+g(h+ηB)\begin{aligned} \text{Bernoulli Potential}: && && p &= \frac{1}{2}(u^2 + v^2) + g (h + \eta_B) = \mathcal{KE} + g(h + \eta_B) \end{aligned}

Vector Identity

uu=(×u)×u+12(uu)\vec{\boldsymbol{u}} \cdot \boldsymbol{\nabla}\vec{\boldsymbol{u}} = (\boldsymbol{\nabla}\times \vec{\boldsymbol{u}})\times \vec{\boldsymbol{u}} + \frac{1}{2}\boldsymbol{\nabla}(\vec{\boldsymbol{u}} \cdot \vec{\boldsymbol{u}})
References
  1. Vallis, G. K. (2019). Essentials of Atmospheric and Oceanic Dynamics. Cambridge University Press. 10.1017/9781107588431
  2. Klöwer, M., Jansen, M. F., Claus, M., Greatbatch, R. J., & Thomsen, S. (2018). Energy budget-based backscatter in a shallow water model of a double gyre basin. Ocean Modelling, 132, 1–11. 10.1016/j.ocemod.2018.09.006
  3. Guillou, F. L., Lahaye, N., Ubelmann, C., Metref, S., Cosme, E., Ponte, A., Sommer, J. L., Blayo, E., & Vidard, A. (2021). Joint Estimation of Balanced Motions and Internal Tides From Future Wide-Swath Altimetry. Journal of Advances in Modeling Earth Systems, 13(12). 10.1029/2021ms002613