# Baroclinic model formulation¶

## Governing equations¶

The three dimensional model solves the Navier-Stokes equations with Boussinesq and hydrostatic assumptions.

The model uses mode-splitting, i.e. the three dimensional horizontal velocity $$\mathbf{u}$$ is split into a depth average $$\bar{\mathbf{u}}$$ and a deviation $$\mathbf{u}' = \mathbf{u} - \bar{\mathbf{u}}$$. We use the 2D shallow water equations (9)- (10) to solve $$\bar{\mathbf{u}}$$, and the 3D momentum equation (4) to solve $$\mathbf{u}'$$. See modules shallowwater_eq and momentum_eq for more information.

Since the model is hydrostatic, the vertical velocity is solved diagnostically from the continuity equation (18). The solver is implemented in VerticalVelocitySolver.

Water temperature and salinity are modeled as tracers by the means of the tracer advection-diffusion equation (13). Options ModelOptions3d.solve_temperature, ModelOptions3d.solve_salinity determine whether the dynamic equations are solved at run time. If not, we treat these variables as constants whose value is set with ModelOptions3d.constant_temperature and ModelOptions3d.constant_salinity options.

In baroclinic simulations the water density, $$\rho$$, depends on the temperature and salinity via the equation of state: $$\rho = \rho'(T, S, p) + \rho_0$$, where $$\rho_0$$ is a constant reference density. The equation of state is set by the option ModelOptions3d.equation_of_state_type, and $$\rho_0$$ is defined in the physical_constants module. Water density affects the internal pressure gradient through the baroclinic head, $$r$$, which we can solve diagnostically from (2). The internal pressure gradient, $$\mathbf{F}_{pg} = g\nabla_h r$$, is computed weakly as separate field. The solver is implemented in InternalPressureGradientCalculator.

Setting option ModelOptions3d.use_baroclinic_formulation=True activates baroclinicity, i.e. the computation of water density, baroclinic head and (both 2D and 3D) internal pressure gradients. If use_baroclinic_formulation=False, water density is not computed. Temperature and salinity may still be simulated, but they are treated as passive tracers.

The following tables summarize the prognostic and diagnostic variables.

Variable

Symbol

Dynamic equation

Thetis field name

Water elevation

$$\eta$$

(9)

elev_2d, elev_3d

Depth av. velocity

$$\bar{\mathbf{u}}$$

(10)

uv_2d

3D velocity

$$\mathbf{u}'$$

(4)

uv_3d

Water temperature

$$T$$

(13)

temp_3d

Water salinity

$$S$$

(13)

salt_3d

Table 1. Prognostic variables in the 3D model.

Variable

Symbol

Equation

Thetis field name

Vertical velocity

$$w$$

(18)

w_3d

Water density

$$\rho$$

(17)

rho_3d

$$r$$

(2)

baroc_head_3d

$$\mathbf{F}_{pg}$$

(3)

int_pg_3d

Table 2. Diagnostic variables in the 3D model.

## Spatial discretization¶

Currently Thetis supports two finite element families: Equal order Discontinuous Galerkin (DG) (option 'dg-dg'), and mimetic Raviart-Thomas-DG family ('rt-dg'). The element family is set by the ModelOptions3d.element_family option. Currently only linear elements are supported, i.e. ModelOptions3d.polynomial_degree must be 1.

The function spaces for both element families are summarized in the following tables.

Variable

Symbol

Function space

Water elevation

$$\eta$$

P1DG

Depth av. velocity

$$\bar{\mathbf{u}}$$

P1DG

3D velocity

$$\mathbf{u}'$$

P1DG x P1DG

Water temperature

$$T$$

P1DG x P1DG

Water salinity

$$S$$

P1DG x P1DG

Vertical velocity

$$w$$

P1DG x P1DG

Water density

$$\rho$$

P1DG x P1DG

$$r$$

P1DG x P2

$$\mathbf{F}_{pg}$$

P1DG x P1DG

Table 3. Equal order Discontinuous Galerkin function spaces (degree=1).

Variable

Symbol

Function space

Water elevation

$$\eta$$

P1DG

Depth av. velocity

$$\bar{\mathbf{u}}$$

RT2

3D velocity

$$\mathbf{u}'$$

HDiv(RT2 x P1DG)

Water temperature

$$T$$

P1DG x P1DG

Water salinity

$$S$$

P1DG x P1DG

Vertical velocity

$$w$$

HDiv(P1DG x P2)

Water density

$$\rho$$

P1DG x P1DG

$$r$$

P1DG x P2

$$\mathbf{F}_{pg}$$

HDiv(RT2 x P1DG)

Table 4. Raviart-Thomas Discontinuous Galerkin function spaces (degree=1).

In both cases the tracers belong to fully discontinuous P1DG x P1DG function space. Tracer advection is solved with upwinding method and slope limiters (see VertexBasedP1DGLimiter).

## Temporal discretization¶

The system of coupled equations is marched in time with a CoupledTimeIntegrator. The time integration method is set by ModelOptions3d.timestepper_type option. Currently supported time integrators are listed below.

Time integrator

Thetis class

2D mode

ALE mesh support

Description

'SSPRK22'

CoupledTwoStageRK

implicit

yes

Coupled method based on SSPRK(2,2) scheme

'LeapFrog'

CoupledLeapFrogAM3

implicit

yes

Table 5. Supported 3D time integrators.

The 2D and 3D time steps can be set via ModelOptions3d.timestep and ModelOptions3d.timestep_2d options. The 2D mode can be treated either implicitly or explicitly. In case of an implicit 2D mode, the 2D time step is equal to the 3D time step and timestep_2d option is ignored.

Thetis can also estimate the maximum stable time step based on the mesh resolution, used element family and time integration scheme. To use this feature, the user should provide the following estimates:

When the simulation initializes, Thetis will compute the maximal feasible time step:

Coupled time integrator: CoupledTwoStageRK
2D time integrator: TwoStageTrapezoid
3D time integrator: SSPRK22ALE
3D implicit time integrator: BackwardEuler
- dt 2d swe: 7.34794172415