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 (17). 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\) (17) w_3d
Water density \(\rho\) (16) rho_3d
Baroclinic head \(r\) (2) baroc_head_3d
Pressure gradient \(\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
Baroclinic head \(r\) P1DG x P2
Pressure gradient \(\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
Baroclinic head \(r\) P1DG x P2
Pressure gradient \(\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
'SSPRK33' CoupledSSPRKSemiImplicit implicit no Coupled method based on SSPRK(3,3) scheme
'LeapFrog' CoupledLeapFrogAM3 implicit yes Leapfrog Adams-Moulton 3 method
'ERKALE' CoupledERKALE explicit yes Fully explicit RK scheme
'IMEXALE' CoupledIMEXALE implicit yes Coupled IMEX method

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
- dt h. advection: 213.200697179
- dt v. advection: 729.166666667
- dt viscosity: 45454.5372777
- CFL adjusted dt: 2D: inf 3D: 213.200697179
- chosen dt: 2D: 213.0 3D: 213.0
- adjusted dt: 2D: 180.0 3D: 180.0