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 (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 advection-diffusion 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 '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 |
|---|---|---|---|---|
|
implicit |
yes |
Coupled method based on SSPRK(2,2) scheme |
|
|
implicit |
yes |
Leapfrog Adams-Moulton 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