# Models¶

Due to computatonnal issues associated with the free surface Navier-Stokes equations, geophysical flows simulations are often carried out with reduced complexity models. With vertically averaged models such as the Saint-Venant system, the water height being an unknown of the system, there is no need to deal with moving meshes. Freshkiss3D’s numerical schemes are based on multilayer Saint-Venant systems with mass-exchanges (i.e. layer-averaged Euler or layer-averaged Navier-Stokes systems) of reduced complexity with respect to the Navier-Stokes equations but enriched with respect to the classic Shallow Water system. Layer-averaged models implemented in Freshkiss3D include topography, viscosity and friction source terms as well as variable density and tracer convection-diffusion.

In this section the Navier-Stokes with free surface system is presented in the introduction and then the focus is put on layer-averaged models implemented in Freshkiss3D. Given the complexity of the numerical methods involved, this documentation only scratches the surface of Freshkiss3D’s solver and it is highly recommended to refer to the References section.

## I - System of equations for free surface flows¶

The fluid domain along \(z\) axis is delimited by free suface \(\eta(t,x,y)\) and topography \(z_b(x,y)\) and water depth is given by \(h(t,x,y)=\eta(t,x,y)-z_b(x,y)\).

### I.2 - Tracer¶

The governing convection-diffusion equation for a tracer \(T\) inside the fluid domain is given by:

Important

**(T) Tracer convection-diffusion equation**:

with \(\nu_T\) is defined by \(D_T/\rho_0\) with \(D_T\) the diffusion coefficient.

Note

In Freshkiss3D

All default parameters are in the

*conf_dict.py*file in the*problem*directory and could be changed (see details in Python API,`problem`

)Tracers are optional and can be added to the simulation with the

`fk.Tracer`

class.By default the density \(\rho\) (see

**Navier-Stokes system**) is set to \(\rho_0\) but can be set as a function of the tracer with the`variable_density`

key of`numerical_parameters`

in the`fk.Problem`

class.The default density state laws included in Freshkiss3D can either take into account the effect of temperature, salinity and particle suspensions. See the parameters of the

`fk.Tracer`

class in the Python API.

### I.3 - Boundary conditions and initial conditions¶

To complete the **Navier-Stokes system with free surface** boundary conditions
and initial conditions have to be provided.
Bottom and free surface kinematic and dynamic boundary conditions are defined in Freshkiss3D
by:

and:

with \({\bf U}_s=U(t,x,y,z=\eta)\), \({\bf U}_b=U(t,x,y,z=Z)\), \({\bf n}_s\) and \({\bf n}_b\) the unit outwards normals at the free surface and bottom, \(p^{atm}\) the atmospheric pressure, \(\tau_w\) the wind shear stress at the free surface and \(\kappa_b(h,\mathbf{u})\) the friction coefficient.

Note

In Freshkiss3D (cf. Python API)

\(\tau_w\) is optional and can be defined with the

`fk.WindEffect`

external effect class.\(\kappa_b(h,\mathbf{u})\) is optional and can be defined with the

`fk.Friction`

class.

On horizontal boundaries several types of boundary condition can be prescribed depending on the flow regime. The way of imposing these conditions depends on the layer-averaged model’s formulation and will not be detailed here. Additionally one has also to prescribe initial conditions to complete the model.

## II - Layer-averaged models¶

Layer-averaged models are the core of Freshkiss3D and allow a 3D modelling of free surface flows. Depending on the initial hypotheses (i.e. user selected options) the models solved by Freshkiss3D behind the curtain can vary. In this section we briefly present some mathematical aspects of the vertical discretization as well as the different models obtained depending on various hypotheses. Some hints on how to select specific options are given with links to the Python API.

### II.1 - Layerwise discretization¶

We consider a discretization of the fluid domain by layers where the layer \(\alpha\) contains the points of coordinates \((x,y,z)\) with \(z \in L_\alpha(t,x,y)= (z_{\alpha-1/2},z_{\alpha+1/2})\) and \(\{z_{\alpha+1/2}\}_{{\alpha} \in \{1,\dots,N \}}\) is defined for \(\alpha \in \{1,\dots,N \}\) by:

Using the previous notations, let us consider the space \(\Bbb{P}_{0}\) of piecewise constant functions defined by

where \(\mathbb{1}_{z \in L_\alpha(t,x,y)}(z)\) is the characteristic function of the layer \(L_\alpha(t,x,y)\). Using this formalism, the projection of \(u\), \(v\), \(w\) and \(T\) on \(\Bbb{P}_{0}\) is a piecewise constant function defined by:

Note that the free surface changing over time, the \(\Bbb{P}_{0}\) space as well as the piecewise constant function projections are time dependent.

### II.2 - The layer-averaged Euler system¶

Let us first neglect the viscosity and assume that the density is constant (i.e. \(\rho =\rho_0\)). Using the space \(\Bbb{P}_0\) and the above decomposition, the Galerkin approximation of the incompressible and hydrostatic Euler equations leads to the system:

Important

**(LE.1) The layer-averaged Euler system**:

for \(\alpha \in \{1,\dots,N \}\):

The quantity \(G_{\alpha+1/2}\) (resp. \(G_{\alpha-1/2}\)) corresponds to the mass exchange accross the interface \(z_{\alpha+1/2}\) (resp. \(z_{\alpha-1/2}\)) and \(G_{\alpha+1/2}\) is defined by:

for \(\alpha=1,\ldots,N\). The velocities at the interfaces \({\bf u}_{\alpha+1/2}\) are defined by:

For smooth solutions, the system **(LE.1)** satisfies an energy balance (see References for complete proof).
Note that the vertical velocity is not an unknown of the system but can be computed considering the
divergence free condition and we obtain the following piecewise approximation of the vertical velocity:

Also note that when the number of layers is equal to one, the classical Saint-Venant system with topography is obtained:

Important

**(SV) The Saint-Venant system**:

Note

In Freshkiss3D (cf. Python API)

**The layer-averaged Euler system is solved by default when no other option is selected**Layers are defined with the

`fk.Layer`

class which is initialized with a user defined number of layers:`NL`

.

Applying the same decomposition to the tracer equation gives:

Important

**(LT) The layer-averaged tracer equation**:

for \(\alpha \in \{1,\dots,N \}\):

where \(G_{\alpha+1/2}\) and \(G_{\alpha-1/2}\) are defined as in **(II.2)**.

### II.3 - The layer-averaged Euler system with variable density¶

As mentioned before, the variable density model relies on the Boussinesq assumption, which means that the variable density only appears in the gravity term. When the density is not constant, integrating the vertical pressure gradient leads to a slighly different model as one cannot assume \(p=\rho_0 g (\eta-z)\) anymore. Now \(p = g \int_{z}^{\eta} \rho dz\). Again, using piecewise approximation one can write the pressure as:

The Galerkin approximation of the incompressible and hydrostatic Euler equations with variable density leads to the system:

Important

**(LE.2) The layer-averaged Euler system with variable density**:

for \(\alpha \in \{1,\dots,N \}\):

where \(G_{\alpha+1/2}\) (resp. \(G_{\alpha-1/2}\)),
\(z_{\alpha+1/2}\) (resp. \(z_{\alpha-1/2}\)),
\({\bf u}_{\alpha+1/2}\) (resp. \({\bf u}_{\alpha-1/2}\)) are defined as in **(LE.1)**,
\(\rho_{\alpha}=f(T_{\alpha})\) with \(f\) a state law and
\(p_{\alpha}\), \(p_{\alpha+1/2}\) are defined by:

Note

In Freshkiss3D (cf. Python API)

As mentionned before, variable density model can be activated by setting the

`variable_density`

key of`numerical_parameters`

to`True`

.The

`ipres`

key of`numerical_parameters`

allows the user to select the pressure term formulation of**(LE.2)**instead of**(LE.1)**even if \(\rho\) is set to constant, in which case \(\rho_{\alpha}\) is replaced in all equations by \(\rho_0\) and the tracer is optional. Note that when`boussinesq`

is prescribed,`ipres`

is automatically set to`True`

and tracer is mandatory.

## III - Numerical scheme¶

The layer-averaged Navier-Stokes system has the form:

where the state vector is defined by:

with \(q_{x,\alpha}=h_\alpha u_{\alpha}\), \(q_{y,\alpha}=h_\alpha v_{\alpha}\). \(F(U)\) stands for the fluxes of the conservative part, and \(S_b(U)\), \(S_e(U)\) and \(S_{v,f}(U)\) are the source terms, representing respectively the topography effects, the momentum exchanges and the viscous and friction effects.

### III.1 - Time discretization¶

Let us denote discrete times by \(t^n\) and \(t^{n+1}=t^n+\Delta t^n\). For the time discretisation of the layer-averaged Navier-Stokes system a time splitting technique is adopted to treat the conservative part, the vertical exchanges and the viscous part successively:

The conservative part is fully explicit and includes the topography source term which is treated with a well-balanced hydrostatic reconstruction scheme. The vertical exchanges terms are semi-implicit, \(G\) being computed with \(U^{n}\). The viscous and friction effects are semi-implicit as well, the horizontal viscosity being fully explicit and the vertical viscosity and friction (treated at the same time) fully implicit. The coriolis force is an semi-implicit.

Noting \(U^{n+1} = U^n + \Delta t^n f(U^n)\) the fully explicit first order scheme, the second order scheme extension implemented in Freshkiss3D, based on a modified version of the Heun scheme, is:

with:

Since \(\gamma \geq 0\), \(U^{n+1}\) is a convex combination of \(U^n\) and \(U^{(2)}\) so the scheme preserves the positivity. For the previous relations \(\Delta t_1^{n}\) and \(\Delta t_2^{n}\) respectively satisfy the CFL conditions associated with \(U^n\) and \(U^{(1)}\).

Note

In Freshkiss3D (cf. Python API)

Time discretization parameters are passed to

`fk.Problem`

class via the`fk.Simutime`

class.By default time step \(\Delta t^n\) is computed with a CFL condition to ensure stability. The default CFL number is set to 0.9 but can be changed with the parameter

`cfl_adjuster`

key of the`numerical_parameters`

dictionary. Although it is not recommended, the user can force a constant time step by passing`constant_dt`

to`fk.Simutime`

.In

`fk.Simutime`

you can set the`second_order`

to`True`

to select second order time scheme.By default semi-implicit methods are used for \(S_e\) and \(S_{v,f}\) but the user can switch back to fully explicit methods by setting

`implicit_exchanges`

and`implicit_vertical_viscosity`

parameters to`False`

in the`fk.Problem`

class.

### III.2 - Space discretization¶

Let \(\Omega\) denote the computational domain with boundary \(\Gamma\), which we assume is polygonal. Let \(T_h\) be a triangulation of \(\Omega\) for which the vertices are denoted by \(P_i\) with \(S_i\) the set of interior nodes and \(G_i\) the set of boundary nodes. For the space discretization of the system, we use an hybrid method. A finite volume technique is used for the conservative part and exchange terms and a \(\Bbb{P}_1\) finite element approach on \(T_h\) for the viscous part.

We recall now the general formalism of finite volumes on unstructured meshes. The dual cells \(C_i\) are obtained by joining the centers of mass of the triangles surrounding each vertex \(P_i\). We use the following notations:

\(K_i\), set of subscripts of nodes \(P_j\) surrounding \(P_i\),

\(|C_i|\), area of \(C_i\),

\(\Gamma_{ij}\), boundary edge between the cells \(C_i\) and \(C_j\),

\(L_{ij}\), length of \(\Gamma_{ij}\),

\({\bf n}_{ij}\), unit normal to \(\Gamma_{ij}\), outward to \(C_i\) (\({\bf n}_{ji}= -{\bf n}_{ij}\)).

If \(P_i\) is a node belonging to the boundary \(\Gamma\), we join the centers of mass of the triangles adjacent to the boundary to the middle of the edge belonging to \(\Gamma\) and we denote

\(\Gamma_i\), the two edges of \(C_i\) belonging to \(\Gamma\),

\(L_i\), length of \(\Gamma_i\) (for sake of simplicity we assume in the following that \(L_i = 0\) if \(P_i\) does not belong to \(\Gamma\)),

\({\bf n}_i\), the unit outward normal defined by averaging the two adjacent normals.

Note

In Freshkiss3D (cf. Python API)

The triangular mesh is passed to the

`fk.Problem`

class via the`fk.TriangularMesh`

class.The

`fk.TriangularMesh`

class can be defined from a file with`fk.TriangularMesh.from_msh_file('file.msh')`

Freshkiss3D mesh IO supports 3 different mesh formats for your mesh file: FreeFem++

`.msh`

, Medit`.mesh`

and Gmsh`.msh`

(format msh2).See the

**Meshes**section of the Gallery of examples to learn how to manipulate meshes with Freshkiss3D

We define the piecewise constant functions \(U^n(x,y)\) on cells \(C_i\) corresponding to time \(t^n\) and \(z_b(x,y)\) as:

with \(U^n_i = (h^n_i,q^n_{x,1,i},\ldots,q^n_{x,N,i},q^n_{y,1,i},\ldots,q^n_{y,N,i})^T\) i.e. \(U^n_i \approx \frac{1}{|C_i|}\int_{C_i} U(t^n,x,y)dxdy\), \(z_{b,i} \approx \frac{1}{|C_i|}\int_{C_i} z_b(x,y)dxdy\). We will also use the notation \(U^n_{\alpha,i} \approx \frac{1}{|C_i|}\int_{C_i}U_\alpha(t^n,x,y)dxdy\). A finite volume scheme for solving the conservative part is:

with \(\sigma_{i,j} = \frac{\Delta t^nL_{i,j}}{|C_i|},\quad \sigma_{i} = \frac{\Delta t^nL_{i}}{|C_i|}\). Here we consider first-order explicit schemes where:

The term \(\cal{F}_{i,j}\) denotes an interpolation of the normal component of the flux \(F(U^n). {\bf n}_{i,j}\) along the edge \(C_{i,j}\). The functions \(F(U^n_i,U^n_j,z^n_{b,i} - z^n_{b,j},{\bf n}_{i,j})\in {\cal R}^{2N+1}\) are the numerical fluxes which are computed using a kinetic interpretation of the system (see References). The source terms discretization (vertical exchanges, viscosity,…) as well as the second order scheme have been omitted but are detailed in reference papers (see References).

Note

In Freshkiss3D (cf. Python API)

In

`fk.Problem`

you can either set`second_order`

parameter or`second_order`

key of`numerical_parameters`

dictionary to`True`

in order to select the second order space scheme.

### III.3 - Initial conditions¶

Initial conditions are needed to define initial state vector:

as well as initial values for topography: \(z^0_b\). By default all variables are set to zero in Freshkiss3D if no initial values are provided.

Note

In Freshkiss3D (cf. Python API):

The topography is initialized when calling the

`fk.Layer`

class presented above.The state vector containing

`U`

,`V`

,`H`

,`HL`

,`QX`

, and`QY`

variables is defined with the`fk.Primitive`

class. Initialization takes place when declaring this class for the first time using`init`

parameters, the types of which can be constant float, vector or function.As mentioned before the tracer

`T`

is defined and initialized in`fk.Tracer`

class.

### III.4 - Boundary conditions¶

In order to complete the numerical scheme, one has to provide boundary conditions, i.e. to define the fluxes \({\cal F}_{e,i} = F(U^n_i,U^n_{e,i},{\bf n}_{i})\).

Note

In Freshkiss3D (cf. Python API), boundary conditions are selected with the following class:

`fk.FluvialFlowRate`

: fluvial regime fluid boundary condition with given flow rate.`fk.FluvialHeight`

: fluvial regime fluid boundary condition with given water depth.`fk.TorrentialOut`

: torrential regime fluid boundary condition with free outlet.`fk.Slide`

: solid wall boundary with slip condition \({\bf u} \cdot {\bf n} = 0\).`fk.DampingOut`

: damp waves with coefficient \(\theta\) ; \(\partial_t \bf X = -\theta ({\bf X} - {\bf X}^0)\).Time varying boundary conditions can be passed to

`fk.FluvialFlowRate`

and`fk.FluvialHeight`

using`fk.TimeDependentFlowRate`

and`fk.TimeDependentHeight`

.