Analytical Solutions for Testing of Baroclinic and Wind-Driven Three-Dimensional Hydrodynamic Models

Cheng He

National Water Research Institute, Environment Canada, Burlington, Ontario, Canada.

DOI: 10.4236/jamp.2024.128171                

Abstract

It is always a challenge for a model developer to verify a three-dimensional hydrodynamic model, especially for the baroclinic term over variable topography, due to a lack of observational data sets or suitable analytical solutions. In this paper, exact solutions for the periodic forcing by surface heat flux and wind stress are given by solving the linearized equations of motion neglecting the rotation, advection and horizontal diffusion terms. The temperature at the bottom is set to a prescribed periodic value and a slip condition on flow is enforced at the bottom. The geometry of the quarter annulus, which has been extensively studied for two- and three-dimensional analytical solutions of unstratified water bodies, is used with a general power law variation of the bottom slope in the radial direction and is constant in the azimuthal direction. The analytical solutions are derived in a cylindrical coordinate system, which describes the three-dimensional fluid field in a Cartesian coordinate system. The results presented in this paper should provide a foundation for studying and verifying the baroclinic term over a varied topography in a three-dimensional numerical model.

Keywords

3D Model, Baroclinic Force, Analytical Solutions, Model Verification

1. Introduction

Various numerical methods for the solution of the Navier-Stokes equations have been applied to the complex issues of flood routing, tidal circulation, storm surges and contaminant transport in coastal areas. Proof of the utility of such methods is often demonstrated by comparison of the computed variables with field observations. However, this type of comparison is often incapable of adequately and accurately verifying the dynamics simulated by the numerical model in the study region [1]. The limitations of this approach are largely due to inadequate data and an incomplete understanding of the behavior of the numerical procedure. Observations of three-dimensional components of flow are rarely available throughout the temporal and spatial domains of interest. Thus, in general, databases are inadequate tools for establishing that an existing numerical model is correctly solving the governing equations. This is particularly true for the vertical velocity component which is usually too small to be measured but is of practical importance to such phenomena as coastal upwelling.

The precision with which a numerical scheme solves the full governing equations should also be established, but because of the nonlinearities in the equations, this is difficult to ascertain precisely. Furthermore, the effect of an irregularly shaped boundary on the accuracy of the numerical solution is generally not completely known, although it is acknowledged to be important.

Due to the lack of suitable analytical solutions for testing the output of three-dimensional numerical models, authors [2]-[5] may often resort to comparing various models with one another. The difficulty with this approach is that it is hard to determine which model is correct. The various sources of error and uncertainty in verification can sometimes be obscured in the numerical solution by the adjustment of parameters such as bathymetry, the eddy coefficient, and the bottom drag coefficients.

Several modelers have developed analytical solutions to aid in the development and verification of various models. Zamani and Bombardelli [6] developed four one-dimensional analytical solutions designed to help code verification; these solutions are able to handle the challenges of the scalar transport equation including nonlinearity and spatiotemporal variability of the velocity and dispersion coefficient, and of the source term. Then, a solution of Burgers’ equation in a novel setup was presented. In order to assist in verifying different terms of one- or two-dimensional Shallow Water equations, Delestre et al. [7] collected a number of analytic solutions scattered throughout the literature of various scientific disciplines. By facing similar challenges in different flow study fields, researchers [8]-[10] developed a wide range of analytical solutions or methods to verify different flow simulations. For an open flow field, it is difficult to properly resolve the boundary condition for an analytical solution. In general, they can be divided into “open basin” solutions with an offshore boundary condition and “closed basin” solutions without an offshore boundary condition. By far the most widely used analytical solutions for the purposes of verification of coastal numerical models are the solutions for quarter-annular geometry. This geometry is particularly useful as a test of a model’s ability to accurately simulate a curved boundary. The analytical solution for the depth-averaged case was developed by Lynch and Gray [11] for constant, linear, and quadratic bathymetry for both steady-state wind setup and periodic tidal response. Their analytical solution was extended to three dimensions and rotation by Lynch and Officer [12]. Lynch and Werner [13] calculated and displayed this analytical solution for tidal oscillations near a circular island where the tidal amplitude at the open ocean boundary varies with position. Muccino et al. [14], who recognized the importance of vertical velocity, gave an expression for the vertical velocity associated with the Lynch and Officer case but with a spatially constant tidal amplitude forcing at the open boundary. Hannah and Wright [15] tested a three-dimensional finite element model with an analytical solution for linear depth-dependent wind-driven flows along a rotating coastline. All of these analytical solutions deal with unstratified water bodies. Once stratification over a sharp variable bathymetry is added, numerical models have considerable difficulty resolving the baroclinic terms, especially in the coordinate system [16], such that analytical solutions are needed in the stratified case as well.

Farrow and Patterson [17] [18] developed analytical solutions for two thermally driven flows for a two-dimensional linearly sloping bottom, neglecting advection but allowing for vertical diffusion. In closed basins, several modelers have exploited wind-driven, flat-bottomed analytical solutions. However, Mass and Lam [19] presented some free modes of oscillation over variable topography, uniform rate of stratification, no friction or diffusion, but linearized advection in a two-dimensional closed basin. Their solutions are interesting as there is a two-way coupling between the mass and velocity fields through advection, unlike the solutions of Farell and Patterson [17] [18] and Fortatuno and Baptista [20] where the density field was prescribed. As forcing and friction are neglected in the solutions of Mass and Lam [19], it is not clear how they can be used in coastal model validation.

In this paper, the analytical solutions were developed aimed at verifying three-dimensional hydrodynamic models, especially for improving the numerical schemes used when calculating the baroclinic terms. The solution starts with solving the vertical temperature diffusion equation in a quarter annulus. This produces a variable density field which, in turn, drives the circulation. Wind forcing is also included in the solution. For simplicity, advection, horizontal diffusion, free surface and rotation terms are neglected. However, heat flux through water surface, bottom and internal friction, wind stress, and variable bathymetry have been incorporated into the equations. A solution for the dynamic steady state with a periodic forcing function is obtained. The vertical velocity is computed from the continuity equation. A general power law dependence of depth with offshore distance allows for the testing of numerical schemes for various degrees of bottom slope. In line with the philosophy that these solutions are useful primarily as tools for model verification rather than solutions to actual field problems, emphasis is placed on periodic solutions and other departures from geophysical realistic flows.

2. Analytical Solution in σ-Coordinate System

A similar quarter annular bathymetry and radial coordinate system as used by Lynch and Officer [11] was adopted in this study and is shown in plan and section views in Figure 1(a) and Figure 1(b), respectively. Because the governing equations, boundary conditions and domain geometry do not depend on the azimuthal direction, variation occurs only in the radial and vertical directions.

Figure 1. A vertical cross-section of the simulation domain grids. In the top panel, the grid is shown in radial coordinates, with the Y-axis representing the water depth and the X-axis representing the radial coordinate. The bottom panel shows the grid in the σ-coordinate after transformation. The Y-axis represents the normalized water depth range from 0 to −1 and the grid becomes uniform as a result.

The stratification is generated by a periodic heat flux applied to the water surface. Here we seek periodic solutions of the form $f\left(r,z,t\right)=Re\left[F_0\left(r,z\right){\text{e}}^{i\omega t}\right]$ for temperature and velocity. The governing equation of temperature in linearized and harmonic form, with the horizontal diffusion and advection terms neglected is:

$i\omega T_0-D\frac{{\partial }^2{T}_0}{\partial z^2}=0$ (1)

at the surface

$D\frac{\partial T_0}{\partial z}=h{F}_0\begin{array}{cccc}& & & \end{array}z=0$ (2)

and at the bottom

$T_0=B_0\begin{array}{cccc}& & & \end{array}z=-h$ (3)

in which $T_0\left(z\right)$ is the complex amplitude of the temperature, $D=N_T{h}^2$ is the diffusivity of temperature, ω is the radian frequency, $h{F}_0$ is the amplitude of the heat flux at the water surface, $h\left(r\right)=h_0{r}^m$ ($h_0$ and m are constants) is the bathymetric relation, r is the radial coordinate, and $B_0$ is the amplitude of the temperature specified at the bottom. $N_T$ is a constant and its reciprocal may be interpreted as the thermal spin-up time.

A spatially varying water depth can be transformed into uniform depth using the σ-coordinate, as shown in Figure 1, which will simplify the numerical formulation.

$\sigma =\frac{z}{h}$ (4)

In a general σ-coordinate system (where σ = 0 at the free surface and $\sigma =-1$ at the bottom), the governing equation of heat diffusion and its boundary conditions become:

$\frac{{\partial }^2{T}_0}{\partial {\sigma }^2}-\frac{i\omega }{N_T}{T}_0=0$ (5)

$N_T{h}^2\frac{\partial T_0}{h\partial \sigma }=h{F}_0\begin{array}{ccccc}& & & & \end{array}\sigma =0$ (6)

$T_0=B_0\begin{array}{ccccc}& & & & \end{array}\sigma =-1$ (7)

The solution of Equation (5) with boundary conditions (6) and (7) is:

$T_0\left(\sigma \right)=M_1\cosh \left(\zeta \sigma \right)-M_2{\text{e}}^{-\zeta \sigma }$ (8)

with

$M_1=\frac{1}{\cosh \zeta }\left(B_0+\frac{F_0}{N_T\zeta }{\text{e}}^{\zeta }\right)$ (9)

$M_2=\frac{F_0}{N_T\zeta }$ (10)

${\zeta }^2=\frac{i\omega }{N_T}$ (11)

Next, the density field can be calculated from the temperature solution. Here, we assume that density is linearly related to temperature as:

$\rho \left(\sigma ,t\right)=a_{T}T\left(\sigma ,t\right)+b$ (12)

where a_T is constant and b is a reference density, so the complex amplitude of the density difference, Δρ, with respect to the reference density at any point, is:

$\Delta \rho \left(\sigma \right)=a_T{T}_0\left(\sigma \right)$ (13)

The complex amplitude of baroclinic pressure, P₀, as a function of depth associated with the density variation is:

$\begin{array}{c}{P}_0=\frac{g}{{\rho }_w}h{\displaystyle \underset{\sigma }{\overset{0}{\int }}\Delta \rho \text{d}\sigma }\\ =\frac{g{a}_T}{{\rho }_w\zeta }h\left[M_2\left(1-{\text{e}}^{-\zeta \sigma }\right)-M_1\sinh \left(\zeta \sigma \right)\right]\end{array}$ (14)

where g is gravity and ${\rho }_w$ is water reference density. The baroclinic term in the radial direction can be computed in σ coordinates:

$\begin{array}{c}\frac{\partial P_0}{\partial r}=\left(\frac{\partial }{\partial r}-\frac{2\sigma }{r}\frac{\partial }{\partial \sigma }\right)P_0\\ =\frac{a_{T}g{h}_0{r}^{m-1}}{{\rho }_w\zeta }\left\{M_2\left[m-\left(m+2\sigma \zeta \right){\text{e}}^{-\zeta \sigma }\right]+M_1[2\sigma \zeta \cosh \left(\zeta \sigma \right)-m\sinh \left(\zeta \sigma \right)\right\}\end{array}$ (15)

After the expression for the baroclinic term has been obtained, the three-dimensional velocity field, driven by baroclinic forcing and wind stress can be computed from the momentum equation. In order to solve the equation analytically, it assumes that the free surface elevation would be zero everywhere at all times. It also assumes that the wind stress has the same frequency as the heat flux through the free surface, and although this solution may not represent any realistic physical flows, it can provide an alternative solution for the purpose of testing numerical schemes and program errors.

The momentum equation, without any advection, rotation, free surface pressure and horizontal diffusion terms is:

$\frac{\partial u}{\partial t}=-\frac{g}{{\rho }_w}\frac{\partial P}{\partial r}+D_v\frac{{\partial }^{2}u}{\partial z^2}$ (16)

where $D_v=N_v{h}^2$ is the vertical eddy viscosity, and $N_v$ is a constant. The horizontal velocity is decomposed by separation of variables with the form $u\left(r,z,t\right)=U_0\left(r,z\right){\text{e}}^{i\omega t}=u_r\left(r\right)u_0\left(z\right){\text{e}}^{i\omega t}$ , where $u_0\left(z\right)$ is the unknown complex amplitude of the horizontal velocity, depending only on the vertical coordinate, z, and $u_r\left(r\right)=r^{m-1}$ . By substituting for u and for the solution to the baroclinic term Equation (15), Equation (16) becomes (in the σ-coordinate system):

$\frac{{\partial }^2{u}_0}{\partial {\sigma }^2}-\frac{i\omega }{N_v}{u}_0=\frac{a_{T}g{h}_0}{{\rho }_w\zeta N_v}\left\{M_2\left[m-\left(m+2\sigma \zeta \right){\text{e}}^{-\zeta \sigma }\right]+M_1\left[2\sigma \zeta \cosh \left(\zeta \sigma \right)-m\sinh \left(\zeta \sigma \right)\right]\right\}$ (17)

with the boundary conditions,

(18)

and

$N_v{h}^2\frac{\partial u_0}{h\partial \sigma }={\tau }_{b}h{u}_0\begin{array}{ccccc}& & & & \end{array}\sigma =-1$ (19)

where $h{\tau }_w$ and $h{\tau }_b$ are the wind stress and the bottom slip coefficients, respectively. The general solution to Equation (17) is:

(20)

where

${\xi }^2=\frac{i\omega }{N_v}$ (21)

and Q is the particular solution which has two forms depending on the Prandtl number ($\xi =\zeta ,\mathrm{Pr}=1$ and $\xi \ne \zeta ,\mathrm{Pr}\ne 1$ ), and A and B are constants to be determined. The details of the particular solutions are given in Appendix A.

In the following calculations, we concentrate on the solution of A and B for case $\xi \ne \zeta$ which represents more general applications. The final result of A and B for the case $\xi =\zeta$ is given in Appendix B.

As an aid to understanding the steps followed in the solution procedure, the expression for $\frac{\partial u_0}{\partial \sigma }$ is given here by differentiating Equation (20) and from Q in Appendix A:

$\begin{array}{c}\frac{\partial u_0}{\partial \sigma }=\xi \left(A{\text{e}}^{\xi \sigma }-B{\text{e}}^{-\xi \sigma }\right)+\frac{G}{\Pi }\{M_1\left\{\left[\left(2-m\right)\zeta -\frac{4{\zeta }^3}{\Pi }\right]\cosh \left(\zeta \sigma \right)+2{\zeta }^2\sigma \sinh \left(\zeta \sigma \right)\right\}\\ -M_2\left[\left(2-m\right)\zeta -2\sigma {\zeta }^2-\frac{4{\zeta }^3}{\Pi }\right]{\text{e}}^{-\zeta \sigma }\}\end{array}$ (22)

where

$\Pi ={\zeta }^2-{\xi }^2$ (23)

$G=\frac{g{a}_T{h}_0}{{\rho }_w\zeta N_v}$ (24)

Substituting (22) into the boundary conditions (18) and (19) yields;

$A\xi -B\xi +G\frac{\left(2-m\right)\zeta -\frac{4{\zeta }^3}{\Pi }}{\Pi }\left(M_1-M_2\right)=\frac{{\tau }_w}{N_v}\begin{array}{ccc}& & \end{array}\sigma =0$ (25)

$\begin{array}{l}A\xi {\text{e}}^{-\xi }-B\xi {\text{e}}^{\xi }+\frac{G}{\Pi }\{M_1\left[2{\zeta }^{2}sh\zeta -\frac{(2-m)\zeta -\frac{4{\zeta }^3}{\Pi }}{\Pi }ch\zeta \right]\\ \begin{array}{c}\\ \stackrel{}{}\\ \end{array}-M_2\left\{\zeta {\text{e}}^{\zeta }\left[\left(2-m\right)+2\zeta -\frac{4{\zeta }^2}{\Pi }\right]\right\}\}\\ =\frac{{\tau }_b}{N_v}\{A{\text{e}}^{-\xi }+B{\text{e}}^{\xi }+G\{\frac{M_1}{\Pi }\left[-2\zeta ch\zeta +\left(\frac{4{\zeta }^2}{\Pi }+m\right)sh\zeta \right]\\ -M_2\left[\frac{m}{{\xi }^2}+\frac{{\text{e}}^{\zeta }}{\Pi }\left(m+\frac{4{\zeta }^2}{\Pi }-2\zeta \right)\right]\}\}\begin{array}{cccc}& & & \end{array}\sigma =-1\end{array}$(26)

From (25) and (26), the constants A and B can be obtained:

$A=B+M_4$ (27)

(28)

in which

$-M_2\left\{\frac{m{\tau }_b}{{\xi }^2{N}_v}+\left(m+\frac{4{\zeta }^2}{\Pi }-2\zeta \right)\frac{{\tau }_b{\text{e}}^{\zeta }}{N_v\Pi }-\left[\left(2-m\right)+2\zeta -\frac{4{\zeta }^2}{\Pi }\right]\frac{\zeta {\text{e}}^{\zeta }}{\Pi }\right\}$ (29)

and

$M_4=\frac{\frac{{\tau }_w}{N_v}-\frac{\left(2-m\right)\zeta -\frac{4{\zeta }^3}{\Pi }}{\Pi }G\left(M_1-M_2\right)}{\xi }$ (30)

To obtain the complete three-dimensional velocity field, the Cartesian vertical velocity, w₀, may be determined from the continuity equation. In three-dimensional cylindrical coordinates, with variation in the azimuthal direction neglected and σ coordinates in the vertical direction, the continuity equation is:

$\frac{-2\sigma h}{r^2}\frac{\partial }{\partial \sigma }\left(r{U}_0\right)+\frac{h}{r}\frac{\partial }{\partial r}\left(r{U}_0\right)+\frac{\partial w_0}{\partial \sigma }=0$ (31)

Integration of Equation (31) vertically over σ from σ to 0 and from σ to −1, and incorporation of the surface boundary conditions, $w_0\left(r,0\right)=0$ and the

bottom boundary conditions $w_0\left(r,-1\right)=-U_0\frac{\partial h}{\partial r}=-m{h}_0{r}^{m-1}{U}_0$ , yields:

$w_0\left(r,\sigma \right)={\displaystyle \underset{0}{\overset{\sigma }{\int }}\left[\frac{2\sigma h}{r^2}\frac{\partial }{\partial \sigma }\left(r{U}_0\right)-\frac{h}{r}\frac{\partial }{\partial r}\left(r{U}_0\right)\right]\text{d}\sigma }+w_0\left(r,0\right)$ (32)

and

$w_0\left(r,\sigma \right)={\displaystyle \underset{-1}{\overset{\sigma }{\int }}\left[\frac{2\sigma h}{r^2}\frac{\partial }{\partial \sigma }\left(r{U}_0\right)-\frac{h}{r}\frac{\partial }{\partial r}\left(r{U}_0\right)\right]\text{d}\sigma }+w_0\left(r,-1\right)$ (33)

Multiplying Equation (32) with (σ + 1) and subtracting Equation (33) multiplied with σ yields:

$\begin{array}{c}{w}_0\left(r,\sigma \right)=\left(\sigma +1\right){R\left(r,\sigma \right)|}_0^{\sigma }-\sigma {R\left(r,\sigma \right)|}_{-1}^{\sigma }+\left(\sigma +1\right)w_0\left(r,0\right)-\sigma w_0\left(r,-1\right)\\ =R\left(r,\sigma \right)-\left(\sigma +1\right)R\left(r,0\right)+\sigma R\left(r,-1\right)+\left(\sigma +1\right)w_0\left(r,0\right)-\sigma w_0\left(r,-1\right)\end{array}$ (34)

In which $R\left(r,\sigma \right)$ , $R\left(r,0\right)$ , $R\left(r,-1\right)$ , $w_0\left(r,0\right)$ , $w_0\left(r,-1\right)$ are:

$\begin{array}{c}R\left(r,\sigma \right)={\displaystyle \int \left[\frac{2\sigma h}{r^2}\frac{\partial }{\partial \sigma }\left(r{U}_0\right)-\frac{h}{r}\frac{\partial }{\partial r}\left(r{U}_0\right)\right]\text{d}\sigma }\\ =h{r}^{m-2}\{A\left(2\xi \sigma -2-m\right)\frac{{\text{e}}^{\xi \sigma }}{\xi }+B\left(2\xi \sigma +2+m\right)\frac{{\text{e}}^{-\xi \sigma }}{\xi }\\ +G\{\frac{M_1}{\Pi }\left[4\zeta {\sigma }^2+\frac{4\left(m+2\right)\zeta }{\Pi }+\frac{{\left(2+m\right)}^2}{\zeta }\right]\cosh \left(\zeta \sigma \right)\\ -\left(\frac{2{\zeta }^2}{\Pi }+1+m\right)4\sigma \sinh \zeta \sigma ]-M_2\{[\frac{4{\zeta }^2}{\Pi }\left(2+2\zeta \sigma +m\right)\\ \begin{array}{c}\\ \end{array}+4\zeta \sigma \left(1+m+\zeta \sigma \right)+{\left(2+m\right)}^2]\frac{{\text{e}}^{-\zeta \sigma }}{\zeta \Pi }-\frac{m^2\sigma }{{\xi }^2}\}\}\}\end{array}$ (35)

$R\left(r,0\right)=h{r}^{m-2}\left(2+m\right)\left\{-\frac{A}{\xi }+\frac{B}{\xi }+G\left\{\frac{M_1}{\Pi }\left[\frac{4\zeta }{\Pi }+\frac{2+m}{\zeta }\right]-\frac{M_2}{\Pi \zeta }\left[\left(2+m\right)+\frac{4{\zeta }^2}{\Pi }\right]\right\}\right\}$ (36)

$\begin{array}{c}R\left(r,-1\right)=h{r}^{m-2}\{-A\left(2\xi +2+m\right)\frac{{\text{e}}^{-\xi }}{\xi }+B\left(-2\xi +2+m\right)\frac{{\text{e}}^{\xi }}{\xi }\\ +G\{\frac{M_1}{\Pi }\left[4\zeta +\frac{{\left(2+m\right)}^2}{\zeta }+\frac{4\zeta \left(m+2\right)}{\Pi }\right]\cosh \zeta -4\left(1+m+\frac{2{\zeta }^2}{\Pi }\right)\sinh \zeta ]\\ -M_2\left\{\left[{\left(2+m\right)}^2-4\zeta \left(1+m-\zeta \right)+\frac{4{\zeta }^2\left(2-2\zeta +m\right)}{\Pi }\right]\frac{{\text{e}}^{\zeta }}{\zeta \Pi }+\frac{m^2}{{\xi }^2}\right\}\}\}\end{array}$ (37)

$w_0\left(r,0\right)=0$ (38)

and

$\begin{array}{c}{w}_0\left(r,-1\right)=-U_0\left(r,-1\right)\frac{\partial h}{\partial r}\\ =-h_{0}m{r}^{m-2}\{A{\text{e}}^{-\xi }+B{\text{e}}^{\xi }+G\{\frac{M_1}{\Pi }\left[-2\zeta \cosh \zeta +\left(\frac{4{\zeta }^2}{\Pi }+m\right)\sinh \zeta \right]\\ -M_2\left[\left(m+\frac{4{\zeta }^2}{\Pi }-2\zeta \right)\frac{{\text{e}}^{\zeta }}{\Pi }+\frac{m}{{\xi }^2}\right]\}\}\end{array}$ (39)

3. Test Cases

The solution presented herein incorporates several features which should be of interest to those working with numerical models. The inclusion of wind stress, baroclinic forcing, vertical mixing, bottom friction, nonlinear bottom slope, and horizontal and vertical velocity gives a broad spectrum of conditions against which model features can be tested. The availability of complete solutions for velocity provides the means for verification of the computed flow field, which is in many practical cases the most important aspect of a problem. The solutions for two-dimensional polar geometry are especially interesting, since they can be used to compare with three-dimensional numerical model outputs in a Cartesian coordinate system.

Here, two analytical evaluations have been given to compare with the output of a three-dimensional finite element numerical model output, in order to show the utility of the analytical solution. At the same time, they also show a sensitivity to vertical viscosity which is of considerable practical importance. The detailed description of the three-dimensional finite element model, LACOM3D, (Lake and Coastal Model in Three-Dimensions), and an extensive comparison with different analytical solutions are left to a future study.

We solve this problem on the mesh depicted in Figure 2(a), with 825 nodes and 1536 triangular elements in the horizontal and 16 evenly spaced levels in the vertical. The bottom slope, as shown in Figure 2(b), is quadratic in r (radius) such that $h=h_0{r}^2$ with $h_0=6.25\times {10}^{-9}/\text{m}$ . At the surface, the heat flux, $F_0=5\times {10}^{-4}\text{C}\cdot {\text{s}}^{-1}$ , and the wind stress, ${\tau }_w=-5\times {10}^{-11}\text{m}\cdot {\text{s}}^{-2}$ , are applied with $\omega =7.27205\times {10}^{-5}{\text{s}}^{-1}$ corresponding to a forcing period of 24 hours. The temperature at the bottom, $B_0=4˚\text{C}$ , and the bottom friction are specified in the examples. The unforced open boundary is applied at both the inner and outer boundaries. The temperature volume expansion coefficient, a_T, is determined by a linear least square fit of the general relationship between density and temperature for freshwater, which is −0.169695 kg·m⁻³ ˚C⁻¹. The diffusivity coefficient, $N_T=1\times {10}^{-5}{\text{s}}^{-1}$ corresponding to a spin-up time of approximately 1 day, is used during the test period. The numerical model was run from rest. The baroclinic term is computed in the σ coordinate system in LACOM3D as opposed to the z-level coordinate system. The time step is 10 seconds.

(a)

(b)

Figure 2. (a) Plan view of the quarter annulus with inner and outer boundaries. (b) Cross-sectional view of the quarter annulus with quadratic bottom dependence.

Figure 3(a) and Figure 3(b) display the results for Pr = 1. Since the strongest wind stress and maximum heat flux at the surface are at t = 0 and at 24-hour intervals, and also because the spin-up time is about 24 hours, in Figure 3(a) the comparison of currents and temperatures along a radial cross-section between the analytical and the numerical model outputs at t = 24 hours is given. It can be seen that after a 24-hour spin-up time, the numerical model output matches the analytical result very well. It also shows that the upper portion is dominated by the wind and the lower portion is driven by the broclinic forcing. Figure 3(b) shows that after 102 hours when the external forces are minimal, the baroclinic term controls most of the water body except near the surface due to inertia. Moreover, after four cycles the numerical model does not show any evidence of error accumulation.

(a)

(b)

Figure 3. (a) Comparison of currents and temperature along a radial cross-section between the analytical and the numerical model output for viscosity equal to temperature diffusivity after 24 hours (when external forces are at their maxima). (b) The same comparison as (a) after 102 hours (as external forces are minimal)

The second test uses the same parameters as the first one except that vertical eddy viscosity is 10 times larger than the temperature diffusivity. The results are shown in Figure 4(a) and Figure 4(b). It is obvious that not only the velocity profile is much weaker, but also they are more uniform due to the higher viscosity. These solutions show sensitivity to vertical viscosity, which is of considerable practical importance, and as a result, they constitute a useful test case.

(a)

(b)

Figure 4. (a) The same comparison as Figure 2(a) except for a viscosity 10 times larger than the temperature diffusivity. (b) The same comparison as Figure 3(a) after 102 hours of the simulated real-time.

4. Conclusions

Analytical solutions for baroclinic and wind forcing with internal and bottom friction and varying depth have been given. The solutions differ depending on whether or not the temperature diffusivity is equal to the vertical viscosity. All terms included in the analytical solution are important factors for any three-dimensional numerical models, so it should be a useful tool for verifying three-dimensional models and evaluating different numerical schemes. Though the analytical solution represents a two-dimensional domain in the radial coordinate system, it does not reduce its usefulness for verifying three-dimensional numerical models as most models are constructed using a Cartesian coordinate system.

An obvious improvement to the analytical solutions obtained in this paper would be to include the free surface term. Moreover, the inclusion of advection and rotation would be valuable. Continual research into the development of analytical solutions for wider ranges of numerical model testing will be pursued in the future.

Acknowledgements

The support received from the National Water Research Institute, Environment Canada is greatly appreciated.

Appendix

Appendix A: Details of the Particular Solution

For ease of presentation, we rewrite Equation (17) as:

$\frac{{\partial }^2{u}_0}{\partial {\sigma }^2}-{\xi }^2{u}_0=\frac{a_{T}g{h}_0}{{\rho }_w\zeta N_v}\left\{M_2\left[m-\left(m+2\sigma \zeta \right){\text{e}}^{-\zeta \sigma }\right]+M_1\left[2\sigma \zeta \cosh \left(\zeta \sigma \right)-m\sinh \left(\zeta \sigma \right)\right]\right\}$ (A1)

The right-hand side can be rearranged in exponential function form as:

$RH=\frac{a_{T}g{h}_0}{{\rho }_w\zeta N_v}\left[m{M}_2-\frac{m{M}_1}{2}{\text{e}}^{\zeta \sigma }+m\left(\frac{M_1}{2}-M_2\right){\text{e}}^{-\zeta \sigma }+M_1\zeta \sigma {\text{e}}^{\zeta \sigma }+\left(M_1-2M_2\right)\zeta \sigma {\text{e}}^{-\zeta \sigma }\right]$ (A2)

The expression (A2) can be classified as three different forms which are:

$RH1=1$ (A3)

$RH2={\text{e}}^{\pm \zeta \sigma }$ (A4)

$RH3=\sigma {\text{e}}^{\pm \zeta \sigma }$ (A5)

The left-hand side of (A1) can be written as:

$LH=\left(D^2-{\xi }^2\right)u_0=f\left(D\right)u_0$ (A6)

where $D=\frac{\partial }{\partial \sigma }$ and $D^2=\frac{{\partial }^2}{\partial {\sigma }^2}$ . The solution for RH1 term is:

$u_1=-\frac{1}{{\xi }^2}$ (A7)

In order to solve RH2 term the following formulas have to be used depending on weather $\xi \ne \zeta$ or $\xi =\zeta$ . For the $\xi \ne \zeta$ case,

$\frac{1}{f\left(D\right)}{\text{e}}^{a\sigma }=\frac{1}{f\left(a\right)}{\text{e}}^{a\sigma }$ (A8)

therefore

$u_2=\frac{1}{D^2-{\xi }^2}{\text{e}}^{\pm \zeta \sigma }=\frac{1}{{\zeta }^2-{\xi }^2}{\text{e}}^{\pm \zeta \sigma }$ (A9)

For the $\xi =\zeta$ case, due to ${\zeta }^2-{\xi }^2=0$ , the formula (A8) is unsuitable. In this case

$\frac{1}{f\left(D\right)}{\text{e}}^{a\sigma }=\frac{\sigma }{f^{\prime }\left(a\right)}{\text{e}}^{a\sigma }$ (A10)

has to be used. So

$u_2=\frac{1}{D^2-{\xi }^2}{\text{e}}^{\pm \zeta \sigma }=\pm \frac{\sigma }{2\zeta }{\text{e}}^{\pm \zeta \sigma }$ (A11)

By applying the formula

$\frac{1}{f\left(D\right)}{\text{e}}^{a\sigma }V\left(\sigma \right)={\text{e}}^{a\sigma }\frac{1}{f\left(D+a\right)}V\left(\sigma \right)$ (A12)

the RH3 term can be solved as:

$\begin{array}{c}{u}_3=\frac{1}{D^2-{\xi }^2}\sigma {\text{e}}^{\pm \zeta \sigma }\\ ={\text{e}}^{\pm \zeta \sigma }\frac{1}{{\left(D\pm \zeta \right)}^2-{\xi }^2}\sigma \\ ={\text{e}}^{\pm \zeta \sigma }\frac{1}{D^2\pm 2D\zeta +{\zeta }^2-{\xi }^2}\sigma \end{array}$ (A13)

In order to solve u₃, the conditions for $\xi \ne \zeta$ and $\xi =\zeta$ have to be considered separately. For the case $\xi \ne \zeta$ , let:

$\frac{1}{D^2\pm 2D\zeta +{\zeta }^2-{\xi }^2}=c_0+c_{1}D+R\left(D^2\right)$ (A14)

$R\left(D^2\right)$ represents the D terms of higher than the first Power. Multiplying two sides of (A14) by $D^2\pm 2D\zeta +{\zeta }^2-{\xi }^2$ and only keeping the D terms of less than the second power, the constants c₀ and c₁ can be computed by comparing the terms of the same power in D on the two sides of (A15),

$1=c_0\left({\zeta }^2-{\xi }^2\right)\pm c_0\zeta 2D+c_1\left({\zeta }^2-{\xi }^2\right)D$ (A15)

$c_0=\frac{1}{{\zeta }^2-{\xi }^2}$ (A16)

$\pm 2c_0\zeta +c_1\left({\zeta }^2-{\xi }^2\right)=0$ (A17)

$c_1=\mp \frac{2\zeta }{{\left({\zeta }^2-{\xi }^2\right)}^2}$ (A18)

So u₃ can be obtained as

$\begin{array}{c}{u}_3={\text{e}}^{\pm \zeta \sigma }\frac{1}{D^2\pm 2D\zeta +{\zeta }^2-{\xi }^2}\sigma \\ ={\text{e}}^{\pm \zeta \sigma }\left(\frac{1}{{\zeta }^2-{\xi }^2}\mp \frac{2\zeta }{{\left({\zeta }^2-{\xi }^2\right)}^2}\frac{\partial }{\partial \sigma }\right)\sigma \\ ={\text{e}}^{\pm \zeta \sigma }\left(\frac{\sigma }{{\zeta }^2-{\xi }^2}\mp \frac{2\zeta }{{\left({\zeta }^2-{\xi }^2\right)}^2}\right)\end{array}$ (A19)

When $\xi =\zeta$ , the above method is unable to resolve c₀ and c₁. However, by rearranging $f\left(D\right)$ in (A13) as:

$f\left(D\right)=D^2\pm 2D\zeta =D\left(D\pm 2\zeta \right)=Dh\left(D\right)$ (A20)

and applying the following formula:

$\frac{1}{f\left(D\right)}V\left(\sigma \right)=\frac{1}{D^k}\frac{1}{h\left(D\right)}V\left(\sigma \right)=\underset{k}{\underbrace{{\displaystyle \int \cdots {\displaystyle \int }}}}\frac{1}{h\left(D\right)}V\left(\sigma \right){\left(\text{d}\sigma \right)}^k$ (A21)

and as before, expand $\frac{1}{D\pm 2\zeta }$ as a polynomial of D: the solution (A13) becomes:

$u_3={\text{e}}^{\pm \zeta \sigma }\frac{1}{D^2\pm 2D\zeta }\sigma$

$\begin{array}{c}={\text{e}}^{\pm \zeta \sigma }\frac{1}{D\left(D\pm 2\zeta \right)}\sigma \\ ={\text{e}}^{\pm \zeta \sigma }{\displaystyle \int \frac{1}{D\pm 2\zeta }\sigma \text{d}\sigma }\\ ={\text{e}}^{\pm \zeta \sigma }{\displaystyle \int \frac{1}{2\zeta }\left(\pm 1-\frac{1}{2{\zeta }^2}\frac{\partial }{\partial \sigma }\right)\sigma \text{d}\sigma }\\ =\frac{{\text{e}}^{\pm \zeta \sigma }\sigma }{4\zeta }\left(\pm \sigma -\frac{1}{\zeta }\right)\end{array}$ (A22)

Combining solutions for RH1, RH2 and RH3 with corresponding factors, the particular solution for Equation (A1) can be obtained. When $\xi \ne \zeta$

$\begin{array}{c}{Q}_{\xi \ne \zeta }=\frac{g{a}_T{h}_0}{{\rho }_w\zeta N_v}[-\frac{m{M}_2}{{\xi }^2}-\frac{m{M}_1{\text{e}}^{\zeta \sigma }}{2\left({\zeta }^2-{\xi }^2\right)}+m\left(\frac{M_1}{2}-M_2\right)\frac{{\text{e}}^{-\zeta \sigma }}{{\zeta }^2-{\xi }^2}\\ +\left(\sigma -\frac{2\zeta }{{\zeta }^2-{\xi }^2}\right)\frac{\zeta M_1{\text{e}}^{\zeta \sigma }}{{\zeta }^2-{\xi }^2}+\left(M_1-2M_2\right)\left(\sigma +\frac{2\zeta }{{\zeta }^2-{\xi }^2}\right)\frac{\zeta {\text{e}}^{-\zeta \sigma }}{{\zeta }^2-\xi {}^2}\\ =\frac{g{a}_T{h}_0}{{\rho }_w\zeta N_v}\{\frac{M_1}{{\zeta }^2-{\xi }^2}\left[2\zeta \sigma \cosh \left(\zeta \sigma \right)-\left(\frac{4{\zeta }^2}{{\zeta }^2-{\xi }^2}+m\right)\sinh \left(\zeta \sigma \right)\right]\\ -M_2\left[\frac{m}{{\xi }^2}+\frac{{\text{e}}^{-\zeta \sigma }}{{\zeta }^2-{\xi }^2}\left(2\zeta \sigma +\frac{4{\zeta }^2}{{\zeta }^2-{\xi }^2}+m\right)\right]\}\end{array}$ (A23)

and when $\xi =\zeta$

$\begin{array}{c}{Q}_{\xi =\zeta }=\frac{g{a}_T{h}_0}{{\rho }_w\zeta N_v}[-\frac{m{M}_2}{{\zeta }^2}-\frac{m{M}_1\sigma {\text{e}}^{\zeta \sigma }}{4\zeta }-m\left(\frac{M_1}{2}-M_2\right)\frac{\sigma {\text{e}}^{-\zeta \sigma }}{2\zeta }\\ +\frac{M_1}{4}\left(\sigma -\frac{1}{\zeta }\right){\text{e}}^{\zeta \sigma }-\left(M_1-2M_2\right)\left(\sigma +\frac{1}{\zeta }\right)\frac{\sigma {\text{e}}^{-\zeta \sigma }}{4}]\\ =\frac{g{a}_T{h}_0}{2{\rho }_w{\zeta }^2{N}_v}\{M_2\left[\left(m+1+\zeta \sigma \right)\sigma {\text{e}}^{-\zeta \sigma }-\frac{2m}{\zeta }\right]\\ \begin{array}{c}\\ \end{array}+M_1\sigma \left[\zeta \sigma \sinh \left(\zeta \sigma \right)-\left(m+1\right)\cosh \left(\zeta \sigma \right)\right]\}\end{array}$ (A24)

Appendix B: Vertical Velocity for the Case $\xi =\zeta$

Due to the fact that the particular solution to Equation (17) is different for the cases $\xi =\zeta$ and $\xi \ne \zeta$ , the constants A and B in solution (20) and the solution for the vertical velocity will be different in the two different cases. Here, the final results for the case $\xi =\zeta$ are given, which can be obtained by a similar method in the case $\xi \ne \zeta$ as mentioned before

$A=B+M_4$ (B1)

$B=\frac{G{M}_3+{\text{e}}^{-\xi }\left(\zeta -\frac{{\tau }_b}{N_v}\right)M_4}{2\left(\xi \sinh \xi +\frac{{\tau }_b}{N_v}\cosh \xi \right)}$ (B2)

$M_4=\frac{\frac{{\tau }_w}{N_{\nu }}-G\left(m+1\right)\left(M_2-M_1\right)}{\xi }$ (B3)

$G=\frac{g{a}_T{h}_0}{2{\rho }_w{\xi }^2{N}_v}$ (B4)

$\begin{array}{c}{M}_3=M_2\left\{\left[\left(m+1\right)-\xi \left(1-m\right)-{\xi }^2+\frac{\left(m+1-\zeta \right){\tau }_b}{N_v}\right]{\text{e}}^{\xi }+\frac{2m{\tau }_b}{N_v\xi }\right\}\\ +M_1\left\{\left[{\xi }^2-\left(m+1\right)-\frac{\left(m+1\right){\tau }_b}{N_v}\right]\cosh \xi +\left[\frac{\zeta {\tau }_b}{N_v}-\left(m-1\right)\zeta \right]\sinh \xi \right\}\end{array}$ (B5)

$\begin{array}{c}{w}_0\left(r,\sigma \right)=\left(\sigma +1\right){R\left(r,\sigma \right)|}_0^{\sigma }-\sigma {R\left(r,\sigma \right)|}_{-1}^{\sigma }+\left(\sigma +1\right)w_0\left(r,0\right)-\sigma w_0\left(r,-1\right)\\ =R\left(r,\sigma \right)-\left(\sigma +1\right)R\left(r,0\right)+\sigma R\left(r,-1\right)+\left(\sigma +1\right)w_0(r,0)-\sigma w_0\left(r,-1\right)\end{array}$ (B6)

$\begin{array}{c}R\left(r,\sigma \right)={\displaystyle \int \left[\frac{2\sigma h}{r^2}\frac{\partial }{\partial \sigma }\left(r{U}_0\right)-\frac{h}{r}\frac{\partial }{\partial r}\left(r{U}_0\right)\right]\text{d}\sigma }\\ =h{r}^{m-2}\{A\left(2\xi \sigma -2-m\right)\frac{{\text{e}}^{\xi \sigma }}{\xi }+B\left(2\xi \sigma +2+m\right)\frac{{\text{e}}^{-\xi \sigma }}{\xi }\begin{array}{c}\\ \\ \\ \end{array}\\ +G\{M_2\left\{\left[\frac{\left(\sigma +\frac{1}{\xi }\right)\left(m+2\right)\left(m+3\right)}{\xi }+{\sigma }^2\left(4+3m+2\xi \sigma \right)\right]{\text{e}}^{-p\sigma }+\frac{2m^2\sigma }{\xi }\right\}\\ +M_1\{\frac{\sigma }{\xi }\left[\left(m+2\right)\left(m+3\right)+2{\xi }^2{\sigma }^2\right]\sinh \left(\xi \sigma \right)\\ \begin{array}{c}\\ \\ \\ \end{array}-\frac{1}{{\xi }^2}\left[\left(4+3m\right){\xi }^2{\sigma }^2+\left(m+2\right)\left(m+3\right)\right]\cosh \left(\xi \sigma \right)\}\}\}\end{array}$ (B7)

$R\left(r,0\right)=h{r}^{m-2}\left[\frac{2+m}{\xi }\left(B-A\right)+\frac{G\left(m+2\right)\left(m+3\right)}{{\xi }^2}\left(M_2-M_1\right)\right]$ (B8)

$\begin{array}{c}R\left(r,-1\right)=h{r}^{m-2}\{-A\left(2\xi +2+m\right)\frac{{\text{e}}^{-\xi }}{\xi }+B\left(-2\xi +2+m\right)\frac{{\text{e}}^{\xi }}{\xi }\begin{array}{c}\\ \\ \\ \end{array}\\ +G\{M_2\left\{\left[\left(4+3m-2\xi \right)+\frac{\left(\frac{1}{\xi }-1\right)\left(m+2\right)\left(m+3\right)}{\xi }\right]{\text{e}}^{\xi }-\frac{2m^2}{\xi }\right\}\\ +M_1\{\frac{1}{\xi }\left[\left(m+2\right)\left(m+3\right)+2{\xi }^2\right]\sinh \xi \end{array}$

$\begin{array}{c}\\ \\ \\ \end{array}-\frac{1}{{\xi }^2}\left[\left(4+3m\right){\xi }^2+\left(m+2\right)\left(m+3\right)\right]\cosh \xi \}\}\}$ (B9)

$w_0\left(r,0\right)=0$ (B10)

$\begin{array}{c}{w}_0\left(r,-1\right)=m{r}^{2\left(m-1\right)}{h}_0\{A{\text{e}}^{-\xi }+B{\text{e}}^{\xi }+G\{M_2\left[\left(\xi -m-1\right){\text{e}}^{\xi }-\frac{2m}{\xi }\right]\\ \begin{array}{c}\\ \end{array}+M_1\left[\left(m+1\right)\cosh \xi -\xi \sinh \xi \right]\}\}\end{array}$ (B11)

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Hackett, B. and Røed, L.P. (1994) Numerical Modeling of the Halten Bank Area: A Validation Study. Tellus A: Dynamic Meteorology and Oceanography, 46, 113-133. https://doi.org/10.3402/tellusa.v46i2.15444
[2]	Tartinville, B., Deleersnijder, E., Lazure, P., Proctor, R., Ruddick, K.G. and Uittenbogaard, R.E. (1998) A Coastal Ocean Model Intercomparison Study for a Three-Dimensional Idealised Test Case. Applied Mathematical Modelling, 22, 165-182. https://doi.org/10.1016/s0307-904x(98)00015-8
[3]	Beletsky, D., O’Connor, W.P., Schwab, D.J. and Dietrich, D.E. (1997) Numerical Simulation of Internal Kelvin Waves and Coastal Upwelling Fronts. Journal of Physical Oceanography, 27, 1197-1215. https://doi.org/10.1175/1520-0485(1997)027<1197:nsoikw>2.0.co;2
[4]	Pinho, J., Ferreira, R., Vieira, L. and Schwanenberg, D. (2014) Comparison between Two Hydrodynamic Models for Flooding Simulations at River Lima Basin. Water Resources Management, 29, 431-444. https://doi.org/10.1007/s11269-014-0878-6
[5]	Zamani, B. and Koch, M. (2020) Comparison between Two Hydrodynamic Models in Simulating Physical Processes of a Reservoir with Complex Morphology: Maroon Reservoir. Water, 12, Article No. 814. https://doi.org/10.3390/w12030814
[6]	Zamani, K. and Bombardelli, F.A. (2013) Analytical Solutions of Nonlinear and Variable-Parameter Transport Equations for Verification of Numerical Solvers. Environmental Fluid Mechanics, 14, 711-742. https://doi.org/10.1007/s10652-013-9325-0
[7]	Delestre, O., Lucas, C., Ksinant, P., Darboux, F., Laguerre, C., Vo, T., et al. (2012) SWASHES: A Compilation of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies. International Journal for Numerical Methods in Fluids, 72, 269-300. https://doi.org/10.1002/fld.3741
[8]	Al-Atawi, N.O. and Mashat, D.S. (2022) A Computational Investigation of the Lid-Driven Cavity Flow. American Journal of Computational Mathematics, 12, 283-296. https://doi.org/10.4236/ajcm.2022.122018
[9]	Haque, S., Matsuo, S. and Kim, H. (2021) Three-Dimensional Structure of Non-Equilibrium Homogeneous Condensation Flow in a Supersonic Rectangular Nozzle. Open Journal of Fluid Dynamics, 11, 55-66. https://doi.org/10.4236/ojfd.2021.112004
[10]	Joshi, Y., Savescu, M., Syed, M. and Blackmore, D. (2021) Interesting Features of Three-Dimensional Discrete Lotka-Volterra Dynamics. Applied Mathematics, 12, 694-722. https://doi.org/10.4236/am.2021.128049
[11]	Lynch, D.R. and Gray, W.G. (1978) Analytic Solutions for Computer Flow Model Testing. Journal of the Hydraulics Division, 104, 1409-1428. https://doi.org/10.1061/jyceaj.0005086
[12]	Lynch, D.R. and Officer, C.B. (1985) Analytic Test Cases for Three‐Dimensional Hydrodynamic Models. International Journal for Numerical Methods in Fluids, 5, 529-543. https://doi.org/10.1002/fld.1650050604
[13]	Lynch, D.R. and Werner, F.E. (1987) Three‐Dimensional Hydrodynamics on Finite Elements. Part I: Linearized Harmonic Model. International Journal for Numerical Methods in Fluids, 7, 871-909. https://doi.org/10.1002/fld.1650070902
[14]	Muccino, J.C., Gray, W.G. and Foreman, M.G.G. (1997) Calculation of Vertical Velocity in Three-Dimensional, Shallow Water Equation, Finite Element Models. International Journal for Numerical Methods in Fluids, 25, 779-802. https://doi.org/10.1002/(sici)1097-0363(19971015)25:7<779::aid-fld586>3.0.co;2-f
[15]	Hannah, C.G. and Wright, D.G. (1995) Depth Dependent Analytical and Numerical Solutions for Wind-Driven Flow in the Coastal Ocean. In: Lynch, D.R. and Davies, A.M., Eds., Quantitative Skill Assessment for Coastal Ocean Models, AGU, 125-152.
[16]	Haney, R.L. (1991) On the Pressure Gradient Force over Steep Topography in Sigma Coordinate Ocean Models. Journal of Physical Oceanography, 21, 610-619. https://doi.org/10.1175/1520-0485(1991)021<0610:otpgfo>2.0.co;2
[17]	Farrow, D.E. and Patterson, J.C. (1993) On the Response of a Reservoir Sidearm to Diurnal Heating and Cooling. Journal of Fluid Mechanics, 246, 143-161. https://doi.org/10.1017/s0022112093000072
[18]	Farrow, D.E. and Patterson, J.C. (1994) The Daytime Circulation and Temperature Structure in a Reservoir Sidearm. International Journal of Heat and Mass Transfer, 37, 1957-1968. https://doi.org/10.1016/0017-9310(94)90335-2
[19]	Maas, L.R.M. and Lam, F.A. (1995) Geometric Focusing of Internal Waves. Journal of Fluid Mechanics, 300, 1-41. https://doi.org/10.1017/s0022112095003582
[20]	Fortunato, A.B. and Baptista, A.M. (1996) Evaluation of Horizontal Gradients in Sigma-Coordinate Shallow Water Models. Atmosphere-Ocean, 34, 489-514.