Baroclinic model formulation¶
Governing equations¶
The three dimensional model solves the NavierStokes equations with Boussinesq and hydrostatic assumptions.
The model uses modesplitting, 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 (11)
(12) to solve \(\bar{\mathbf{u}}\), and the
3D momentum equation (5) 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 (21).
The solver is implemented in VerticalVelocitySolver
.
Water temperature and salinity are modeled as tracers by the means of the
tracer advectiondiffusion equation (15).
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 (3).
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\) 


Depth av. velocity 
\(\bar{\mathbf{u}}\) 


3D velocity 
\(\mathbf{u}'\) 


Water temperature 
\(T\) 


Water salinity 
\(S\) 

Table 1. Prognostic variables in the 3D model.
Variable 
Symbol 
Equation 
Thetis field name 

Vertical velocity 
\(w\) 


Water density 
\(\rho\) 


Baroclinic head 
\(r\) 


Pressure gradient 
\(\mathbf{F}_{pg}\) 

Table 2. Diagnostic variables in the 3D model.
Spatial discretization¶
Currently Thetis supports two finite element families:
Equal order Discontinuous Galerkin (DG)
(option 'dgdg'
), and mimetic RaviartThomasDG family ('rtdg'
).
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. RaviartThomas 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 


implicit 
yes 
Coupled method based on SSPRK(2,2) scheme 


implicit 
yes 
Leapfrog AdamsMoulton 3 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:
ModelOptions3d.
horizontal_velocity_scale
: Maximal horizontal velocity scaleModelOptions3d.
vertical_velocity_scale
: Maximal vertical velocity scaleModelOptions3d.
horizontal_viscosity_scale
: Maximal horizontal viscosity scale
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