2. Theoretical framework

Under the Boussinesq approximation, conservation of mass and the Reynolds-averaged Navier–Stokes equations are

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}

(2.2)\begin{gather}\frac{\partial u}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}){u}-f v={-}\frac{1}{\rho_0} \frac{\partial P}{\partial x}+ \boldsymbol{\nabla}\boldsymbol{\cdot}(\nu_h \boldsymbol{\nabla}_h u)+\frac{\partial}{\partial z}\left(\nu_v \frac{\partial u}{\partial z}\right), \end{gather}

(2.3)\begin{gather}\frac{\partial v}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}){v}+f u={-}\frac{1}{\rho_0} \frac{\partial P}{\partial y}+ \boldsymbol{\nabla}\boldsymbol{\cdot}(\nu_h \boldsymbol{\nabla}_h v)+\frac{\partial}{\partial z}\left(\nu_v \frac{\partial v}{\partial z}\right), \end{gather}

(2.4)\begin{gather}\frac{\partial w}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}){w}={-}\frac{1}{\rho_0} \frac{\partial P}{\partial z}+ \boldsymbol{\nabla}\boldsymbol{\cdot}(\nu_h \boldsymbol{\nabla}_h w)+\frac{\partial}{\partial z}\left(\nu_v \frac{\partial w}{\partial z}\right)- \frac{\rho}{\rho_0}g, \end{gather}

where $u,v,w,P,f$, $\nu _h$ and $\nu _v$ are velocity in x, y and z directions, pressure, the Coriolis frequency, and the horizontal and vertical eddy viscosities, respectively. Pressure and density are partitioned into reference, background and perturbation components, and velocities are decomposed into BT and BC components

(2.5)\begin{gather} \rho=\rho_0+\rho_b+\rho', \end{gather}

(2.6)\begin{gather}P=P_0+P_b+P'=g \int_{z}^{\eta} \left( \rho_0+\rho_b(z')+\rho'(z')\right) {\rm d}z', \end{gather}

(2.7)\begin{gather}\boldsymbol{U_h}=\frac{1}{H}\int_{{-}d}^{\eta} \boldsymbol{u_h}= \frac{1}{H} \overline{\boldsymbol{u_h}}\,{\rm d}, \end{gather}

(2.8)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}=0, \end{gather}

(2.9)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u'}=0, \end{gather}

where the $\boldsymbol {h}$ subscript denotes the horizontal ($x,y$ and $u,v$), $H=\eta +d$ is the water column depth ($\eta$ and $-d$ are sea surface height and ocean bottom, respectively), the bar operator provides the depth-integrated value of a parameter ($\bar {\psi }=\int _{-d}^{\eta }\psi \,{\rm d} z$) and capital letter velocities ($U$,$V$) are the BT components, while the prime velocities ($u'$,$v'$) are the BC constituents. The BT and BC vertical velocities can be found using the continuity equation $W=-\boldsymbol {\nabla }_h\boldsymbol {\cdot }(\boldsymbol {U}_h (d+z))$ and $w'=w-W$, respectively.

By using the momentum equations and applying boundary conditions (full derivation is provided in Kang (Reference Kang2010) and summarized in Appendix A), we get the depth-integrated BC kinetic energy ($E'_K$) budget

(2.10)\begin{align} \frac{\partial \overline{E'_K}}{\partial t}+\boldsymbol{\nabla_h} \boldsymbol{\cdot}(\overline{\boldsymbol{u_h}E'_K})+\boldsymbol{\nabla_h} \boldsymbol{\cdot}(\overline{\boldsymbol{u'_h}E_{K_{BT-BC}}}) &={-}\boldsymbol{\nabla_h} \boldsymbol{\cdot}(\overline{\boldsymbol{u_h'}P'})-\overline{\rho'gw'}+A_h \nonumber\\ &\quad +\overline{\boldsymbol{\nabla_h}\boldsymbol{\cdot}(\nu_H\boldsymbol{\nabla_h}E'_K)} -D'-\overline{\epsilon_K}. \end{align}

where $A_h, D', \epsilon_K$ are unclosed, drag and dissipation terms respectively.

For an incompressible density-stratified fluid, the available potential energy ($E_A$) is equal to the difference between the potential energy of the perturbed water column ($E_P$) and the minimum energy attainable thorough adiabatic motion ($E_B$) (Lorenz Reference Lorenz1955; Winters et al. Reference Winters, Lombard, Riley and D'Asaro1995). The concept of $E_A$ has been widely used to study energy budgets in both deep and shallow water systems (Venayagamoorthy & Fringer Reference Venayagamoorthy and Fringer2005; Klymak et al. Reference Klymak, Moum, Nash, Kunze, Girton, Carter, Lee, Sanford and Gregg2006; Scotti, Beardsley & Butman Reference Scotti, Beardsley and Butman2006; Kang & Fringer Reference Kang and Fringer2010, Reference Kang and Fringer2012). For any arbitrary density stratification, $E_A$ can be found from (Kang & Fringer Reference Kang and Fringer2010, Reference Kang and Fringer2012)

(2.11)\begin{equation} E_A=E_P-E_B=g\int_{z-\delta}^{z} \left(\rho(z)-\rho_b(z')-\rho_0 \right) {\rm d}{z'}, \end{equation}

where $\delta$ is the BT and BC vertical displacement of isopycnals. For linear IWs in linear density-stratified media, $E_A$ can further be approximated as (Gill Reference Gill1982; Scotti et al. Reference Scotti, Beardsley and Butman2006; Lamb Reference Lamb2007)

(2.12)\begin{equation} E_A=\frac{\rho_0}{2}N^2\delta^2=\frac{g^2 \rho'^2}{2\rho_0 N^2}, \end{equation}

where $N^2=-({g}/{\rho _0})({\partial \rho }/{\partial z})$ is the squared Brunt–Väisälä frequency. The $E_A$ inside a system can be dissipated, advected or converted to kinetic energy (disregarding the spatio-temporal background density variation ${\partial \rho _b}/{\partial x}, {\partial \rho _b}/{\partial y}, {\partial \rho _b}/{\partial t}\approx 0$)

(2.13)\begin{equation} \overline{\frac{\partial E_A}{\partial t}}+\overline{\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}E_A}=\overline{g w \rho'(z)}+\overline{\epsilon_A}. \end{equation}

where $\epsilon_A$ is the diffusion of potential energy.

By adding (2.10) and (2.13) and tidally averaging, we get

(2.14)\begin{align} & \boldsymbol{\nabla_h}\boldsymbol{\cdot}\left\langle \overline{\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}E_A}+ \overline{\boldsymbol{u_h}E'_k}+\overline{\boldsymbol{u'_h}EK_{BT-BC}}+ \overline{\boldsymbol{u_h'}P'}-\overline{\nu_H\boldsymbol{\nabla_h}E'_K}\right\rangle= \langle\overline{\rho'gw}-\overline{\rho'gw'}\rangle \nonumber\\ &\quad +\langle A_h\rangle-\langle D'\rangle\langle\overline{\epsilon_A-\epsilon_K}\rangle, \end{align}

where the angle bracket operator provides the tidally averaged value of a parameter over a period of T ($\langle \psi \rangle =({1}/{T})\int _{0}^{T} \psi \,{\rm d} t$). The energy conversion rate from BT to BC modes can be found by the summation of $\rho 'gw$ and $-\rho 'gw'$ (Kang & Fringer Reference Kang and Fringer2012). In such a formulation, $C$ consists of both topographic and residual conversion. BTH of isopycnals alone (without the presence of a sloping bottom) does not generate IWs nor convert energy from BT to BC; yet, vertical acceleration over a sloping topography causes a phase lag between the density perturbation and the vertical velocity, $\phi _{\rho ',W}$ (where $\phi _{a,b}$ denotes the phase difference between a and b). The conversion that occurs due to the phase lag between the density perturbation and the vertical velocity can be understood by considering two idealized sinusoidal functions assigned to $\rho '$ and $W$ with period, amplitude and initial phase of $\omega$, $a$ and $b$

(2.15)\begin{equation} \langle C \rangle=\langle \rho'gW \rangle=g\left\langle a_i\sin(\omega t+b_i) a_j\sin(\omega t+b_j)\right\rangle=g\frac{a_ia_j}{2}\cos(\phi_{b_i,b_j}) \end{equation}

leading to positive (negative) conversion rates for phase differences $<90^\circ$ ($>90^\circ$) (Zilberman et al. Reference Zilberman, Becker, Merrifield and Carter2009).

To better understand how BTH contributes to conversion, we consider a system governed by tidal flows with a flat bathymetry in which isopycnals heave up (down) during flood (ebb) tides and $\phi _{W, \rho '} = 90^\circ$. Figure 1 depicts the characteristics of such a system where there are non-zero density perturbations (red curve in figure 1c) and BT vertical velocities (figure 1d) through time. The density perturbations captured by CBD methods is completely induced by BTH (compare panels (b) and (c) in figure 1). Due to the failure of BTH filtering in $\rho _B$ a non-zero $C$ emerges with a zero tidally averaged value (red curve in figure 1f) which is fully BT. By applying a proper time-varying $\rho _b$, the BT displacement of the background density is removed (blue curve in figure 1c); therefore, the residual conversion is eliminated (blue curve in figure 1f).

Figure 1.

Times series of tidal flow over flat bottom showing (a) tidal height, (b) background density, (c) density perturbation, (d) BT vertical velocity at mid-depth, (e) BT conversion rate where $\varDelta \rho _b$ is difference between TVBD and CBD background density and (f) BT–BC conversion rate.

In CBD methods (e.g. both LSY and KF), $\rho _b$ is only a function of $z$; while in TVBD, $\rho _b$ is a function of $z$ and $t$. The temporal dependency of $\rho _b$ creates a new cross-term in the $E_A$ budget which can remove arbitrary vertical displacements (a similar methodology was used by Lu et al. Reference Lu, Wright and Brickman2001; Kunze et al. Reference Kunze, Rosenfeld, Carter and Gregg2002; Zilberman et al. Reference Zilberman, Merrifield, Carter, Luther, Levine and Boyd2011). By removing the BT displacement at any depth we can remove the effects of BTH. To do so, we assume that BT tidal height varies linearly over the water column as $\gamma (x,y,z,t)=\eta (1-{z}/{H})$.

The contribution of nonlinear advection terms in the $E_A$ budget in the formation of $\rho ' g W$ is acknowledged by Kang & Fringer (Reference Kang and Fringer2012). Accounting for nonlinear advection terms results in an additional term in the $E_A$ budget (see Appendix A for full derivation)

(2.16)\begin{equation} \overline{\frac{\partial E_A}{\partial t}}+\overline{\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}E'_K} = \overline{\rho'g w }-\overline{g \int_{z-\delta}^{z}\frac{\partial \rho_b(z',t)}{\partial t} {\rm d} z'}+ \overline{\epsilon_A}. \end{equation}

The total energy budget (kinetic+potential), where $\langle \overline {\rho ' gW-g \int _{z-\delta }^{z}({\partial \rho _b(z',t)}/{\partial t} \,{\rm d} z'})\rangle$ is the modified conversion, then becomes

(2.17)\begin{align} & \boldsymbol{\nabla_h}\boldsymbol{\cdot}\left\langle \overline{\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}E_A}+\overline{\boldsymbol{u_h}E'_k}+ \overline{\boldsymbol{u'_h}EK_{BT-BC}}+\overline{\boldsymbol{u_h'}P'}- \overline{\nu_H\boldsymbol{\nabla_h}E'_K}\right\rangle \nonumber\\ &\quad =\left\langle\overline{\rho'gW-g \int_{z-\delta}^{z} \frac{\partial \rho_b(z',t)}{\partial t} {\rm d} z'}\right\rangle +\langle A_h\rangle-\langle D'\rangle+\langle\overline{\epsilon_A-\epsilon_K}\rangle. \end{align}

This formulation can be considered a modified version of Kang & Fringer (Reference Kang and Fringer2012). The contribution of the non-hydrostatic pressure ($-\overline {({\partial q}/{\partial z})W}$ where $q$ is the non-hydrostatic pressure constituent) can be added directly to this formulation (Kang & Fringer Reference Kang and Fringer2012). The new term in (2.17) can isolate the conversion due to the interaction of IWs and BTH ($IW \times BTH$) and even between IWs ($IW \times IW$) if the background density is formulated accordingly.

3. Numerical set-up

To test the TVBD method and the resulting effects on $C$, we simulated IW generation using SUNTANS (Stanford unstructured non-hydrostatic terrain following averaged Navier–Stokes simulator). SUNTANS is a finite volume model that solves 3-D non-hydrostatic, nonlinear Navier–Stokes equations on unstructured–triangular horizontal grids with $z$-level vertical layers (Fringer, Gerritsen & Street Reference Fringer, Gerritsen and Street2006). SUNTANS has been widely used to study non-hydrostatic, nonlinear phenomena such as IWs and circulation in shallow coastal systems and estuaries (Jachec et al. Reference Jachec, Fringer, Street and Gerritsen2007; Zhang, Fringer & Ramp Reference Zhang, Fringer and Ramp2011; Kang & Fringer Reference Kang and Fringer2012; Nelko, Saha & Chua Reference Nelko, Saha and Chua2014; Xu & Chua Reference Xu and Chua2016). We begin our analysis with two idealized 2-D cases with high spatio-temporal resolution to highlight the TVBD method and compare it to other common methods. To validate the formulation in a more realistic set-up, we compared the results for an idealized 3-D ridge. The bathymetry for the 2-D and the 3-D set-ups can be seen in figures 2 and 3, respectively.

Figure 2.

Vertical velocity at time ${t}/{T}=0.4$ for (a) shallow (case 1) and (b) deep (case 2) scenarios. The generation site over the topographic feature (green circle) and a site of IW reflection from the ocean surface (purple circle) are shown.

Figure 3.

Ideal ridge bathymetry map (a, case 3) and its transects at (b) $X= 300$ km and (c) $Y= 50$ km.

For all cases, we applied a constant background horizontal turbulent viscosity of 1 m$^{2}$ s$^{-1}$ and employed a Mellor–Yamada 2.5 (MY 2.5) turbulence closure scheme for the vertical turbulent viscosity (Mellor & Yamada Reference Mellor and Yamada1982). The horizontal viscosity was set to the minimum value needed to remove numerical oscillations between grid points (10 m resolution). The shortest wavelengths possible in the domain are $O$(60 m), so this did not affect IW generation. We also ran cases with $\nu _h = 0.1$, 0.01 and 0.001 m$^2$ s$^{-1}$ and constant $\nu _v = 0.001$ and 0.0001 m$^2$ s$^{-1}$. Results of all cases were consistent, we therefore only report cases with $\nu_h = 1$ m$^{2}$ s$^{-1}$ and the MY 2.5 closure scheme.

To analyse the effect of depth (shallow vs deep) on $C$, we used idealized functions for density stratification (figure 4) and bathymetry

(3.1)\begin{equation} d=H_0 - 0.5(H_0-h_0)\left(1+\tanh\left(\frac{x-X_{mid}}{L_s}\right)\right), \end{equation}

where $h_0$ and $H_0$ are constants (minimum and maximum depths in the domain), $d$ is the depth at any specific point, $X_{mid}$ and $L_s$ are bathymetric shape constants.

Figure 4.

Background density and Brunt–Väisälä frequency at the beginning of the last 3 tidal cycles for (a) shallow (case 1) and (b) deep (cases 2 and 3) cases.

For all cases, $M_2$ tidal currents ($u=U_0 \sin (\omega t)$) were applied at the ocean-side boundary barotropically and the shore-side boundary was closed. To avoid IW reflection from offshore and onshore boundaries, $u'$ was relaxed using a sponge layer that extended 20 km into the domain at each side. The Coriolis frequency was set to $0.8 \times 10^{-4}$ s$^{-1}$ representing ${\sim }33 ^{\circ }$N latitude. We chose this mid-latitude as it represents a value close to the region (central California) that motivated this study. However, our results were again consistent regardless of the value of $f$ used based on additional simulations.

We examined three case scenarios (shallow 2-D shelf, deep 2-D shelf, 3-D ridge) with different physical inputs (tidal current velocity, bathymetry, eddy viscosity) and environmental parameters (domain size and spatio-temporal resolution) as summarized in table 3 in Appendix A. The models were run for 15 tidal cycles and all analyses were performed on the last 3 tidal cycles. Model set-up files and analysis scripts are available on GitHub at https://github.com/somidvar/suntans/tree/master/TVBDPaperSourceCode.