\(k-\varepsilon\)

\(k-\varepsilon\) closure parameters, states and computation functions.

This module contains the implementation of the \(k-\varepsilon\) model described as a GLS case as in [1] as a Closure instance. The model was traduced from Frotran to JAX with the work of Florian Lemarié and Manolis Perrot [2], the translation was done in part using the work of Anthony Zhou, Linnia Hawkins and Pierre Gentine [3]. The parameters of the closure are available in the KepsParameters class, the closure state in KepsState class. The function keps_step compute one time-step of the closure, which means that it computes the eddy-diffusivity and viscosity. The module contains other functions that are used by this main one. These classes and the function step can be obtained by the prefix tunax.closures.k_epsilon or directly by tunax.closures.

References

[1]

L. Umlauf and H. Burchard. A generic length-scale equation for geophysical turbulence models (2003). Journal of Marine Research 61 pp. 235-265. doi : 10.1357/002224003322005087.

[2]

M. Perrot and F. Lemarié. Energetically consistent Eddy-Diffusivity Mass-Flux convective schemes. Part I: Theory and Models (2024). url : hal.science/hal-04439113.

[3]

A. Zhou, L. Hawkins and P. Gentine. Proof-of-concept: Using ChatGPT to Translate and Modernize an Earth System Model from Fortran to Python/JAX (2024). url : arxiv.org/abs/2405.00018.

class KepsParameters

Parameters and constants for \(k-\varepsilon\).

The first 19 attributes are the parameters of \(k-\varepsilon\) that may be calibrated. This class also contains some physical constants used in the closure computing. The next 16 attirbutes are some physical parameters for k-epsilon that can be modified but that are not specially supposed to be modified. The last 19 attributes are the one for the stability function that are computed from the parameters of \(k-\varepsilon\). The constructor of the class takes as parameters the 19 parameters of the closure and the 16 physical parameters, but not the last 19 stability functions for parameters.

c1

\(k-\varepsilon\) parameter \(c_1\) for the dissipation of the corelation tensor pressure/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=5.

c2

\(k-\varepsilon\) parameter \(c_2\) for the dissipation of the corelation tensor pressure/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=0.8

c3

\(k-\varepsilon\) parameter \(c_3\) for the dissipation of the corelation tensor pressure/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=1.968

c4

\(k-\varepsilon\) parameter \(c_4\) for the dissipation of the corelation tensor pressure/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=1.136

c5

\(k-\varepsilon\) parameter \(c_5\) for the dissipation of the corelation tensor pressure/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=0.

c6

\(k-\varepsilon\) parameter \(c_6\) for the dissipation of the corelation tensor pressure/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=0.4

cb1

\(k-\varepsilon\) parameter \(c_{b1}\) for the dissipation of the corelation tensor buoyancy/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=5.95

cb2

\(k-\varepsilon\) parameter \(c_{b2}\) for the dissipation of the corelation tensor buoyancy/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=.6

cb3

\(k-\varepsilon\) parameter \(c_{b3}\) for the dissipation of the corelation tensor buoyancy/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=1.

cb4

\(k-\varepsilon\) parameter \(c_{b4}\) for the dissipation of the corelation tensor buoyancy/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=0.

cb5

\(k-\varepsilon\) parameter \(c_{b5}\) for the dissipation of the corelation tensor buoyancy/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=0.3333

cbb

\(k-\varepsilon\) parameter \(c_{bb}\) for the dissipation of the corelation tensor buoyancy/velocity (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=.72

c_mu0

\(k-\varepsilon\) parameter \(c_\mu^0\) which links the mixing length to the dissipation (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=0.5477

sig_k

\(k-\varepsilon\) parameter \(\sigma_k\) Schmit number for the dissipation of TKE (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=1.

sig_eps

\(k-\varepsilon\) parameter \(\sigma_\varepsilon\) Schmit number for the dissipation of \(\varepsilon\) (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=1.3

c_eps1

\(k-\varepsilon\) parameter \(c_{\varepsilon 1}\) correction of the \(\varepsilon\) equation (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=1.44

c_eps2

\(k-\varepsilon\) parameter \(c_{\varepsilon 2}\) correction of the \(\varepsilon\) equation (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=1.92

c_eps3m

\(k-\varepsilon\) parameter \(c_{\varepsilon 3}^-\) correction of the \(\varepsilon\) equation (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=-0.4

c_eps3p

\(k-\varepsilon\) parameter \(c_{\varepsilon 3}^+\) correction of the \(\varepsilon\) equation (Umlauf and Burchard notations) \([\text{dimensionless}]\).

Type:

float, default=1.

chk_grav

Charnock coefficient times gravity \([\text{dimensionless}]\).

Type:

float, default=1400.

galp

Parameter for Galperin mixing length limitation \([\text{dimensionless}]\).

Type:

float, default=0.53

z0s_min

Minimal surface roughness length \([\text m]\).

Type:

float, default=1e-2

z0b_min

Minimal bottom roughness length \([\text m]\).

Type:

float, default=1e-4

z0b

Bottom roughness length \([\text m]\).

Type:

float, default=1e-14

akt_min

Minimal and initialization value of eddy-diffusivity \(\left[\text m^2 \cdot \text s^{-1} \right]\).

Type:

float, default=1e-5

akv_min

Minimal and initialization value of eddy-viscosity \(\left[\text m^2 \cdot \text s^{-1} \right]\).

Type:

float, default=1e-4

tke_min

Minimal and initialization value of turbulent kinetic energy (TKE) \(\left[\text m^3 \cdot \text s^{-2} \right]\).

Type:

float, default=1e-6

eps_min

Minimal and initialization value of TKE dissipation \(\left[\text m^2 \cdot \text s^{-3} \right]\).

Type:

float, default=1e-12

c_mu_min

Minimal and initialization value of \(c_\mu\) in GLS formalisim \([\text{dimensionless}]\).

Type:

float, default=0.1

c_mu_prim_min

Minimal and initialization value of c_mu’ in GLS formalisim \([\text{dimensionless}]\).

Type:

float, default=0.1

dir_sfc

Apply a Dirichlet boundary condition at the surface for TKE, else apply a Neumann boundary condition.

Type:

bool, default=False

dir_btm

Apply a Dirichlet boundary condition at the bottom for TKE, else apply a Neumann boundary condition.

Type:

bool, default=True

gls_p

GLS coefficient \(p\) to define \(k-\varepsilon\) \([\text{dimensionless}]\).

Type:

float, default=3

gls_m

GLS coefficient \(m\) to define \(k-\varepsilon\) \([\text{dimensionless}]\).

Type:

float, default=1.5

gls_n

GLS coefficient \(n\) to define \(k-\varepsilon\) \([\text{dimensionless}]\).

Type:

float, default=-1

sf_d0

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_d1

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_d2

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_d3

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_d4

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_d5

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_n0

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_n1

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_n2

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_nb0

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_nb1

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

sf_nb2

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

lim_am0

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

lim_am1

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

lim_am2

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

lim_am3

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

lim_am4

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

lim_am5

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

lim_am6

Limitation coefficient for \(k-\varepsilon\) computed from the parameters \([\text{dimensionless}]\).

Type:

float (not a parameter, computed from the above attributes)

class KepsState

Define the state of the water column for the \(k-\varepsilon\) model.

The first initilisation is done from the minimal values of the different variables given in an instance of :class:’KepsParameters`.

Parameters:

• grid (Grid) – cf. attribute.

• keps_params (KepsParameters) – Used to define the initialization values of the variables.

grid

Geometry of the water column, should be the same than for the State instance used in the model.

Type:

Grid

akt

Eddy-diffusivity on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-1}\right]\).

Type:

float Array of shape (nz+1)

akv

Eddy-viscosity on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-1}\right]\).

Type:

float Array of shape (nz+1)

tke

Turbulent kinetic energy (TKE) denoted \(k\) on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-2}\right]\).

Type:

float Array of shape (nz+1)

eps

TKE dissipation denoted \(\varepsilon\) on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-3}\right]\).

Type:

float Array of shape (nz+1)

c_mu

\(c_\mu\) in GLS formalisim on the interfaces of the cells \([\text{dimensionless}]\).

Type:

float Array of shape (nz+1)

c_mu_prim

\(c_\mu'\) in GLS formalisim on the interfaces of the cells \([\text{dimensionless}]\).

Type:

float Array of shape (nz+1)

keps_step(state, keps_state, dt, keps_params, case_tracable)

Run one time-step of the \(k-\varepsilon\) closure.

The purpose of this function is to get the eddy-diffusivity and eddy-viscosity at the next time- step. It works in 3 steps

1. The Brunt–Väisälä frequency and the shear is computed from the state and the boundary conditions are computed.

2. The equations on \(k\) and \(\varepsilon\) are solved and their values are computed for the next time step.

3. The eddy-diffusivity and viscosity are computed as diagnostic variables and the \(keps_state\) is updated.

Parameters:

• state (State) – Current state of the water column.

• keps_state (KepsState) – Current state of the water column for the variables used by \(k-\varepsilon\).

• dt (float) – Time-step of the forward model \([\text s]\).

• keps_params (KepsParameters) – Values of the parameters used by \(k-\varepsilon\) (time-independant).

• case_tracable (CaseTracable) – Physical parameters and forcings of the model run.

Returns:

keps_state – State of the water column for the variables used by \(k-\varepsilon\) at the next time-step.

Return type:

KepsState

compute_eos(state, case)

Compute density anomaly and Brunt–Väisälä frequency.

Prognostic computation via linear Equation Of State (EOS) :

\(\rho = \rho_0(1-\alpha (T-T_0) + \beta (S-S_0))\)

\(N^2 = - \dfrac g {\rho_0} \partial_z \rho\)

Parameters:

• state (State) – Current state of the water column.

• case (CaseTracable) – Physical parameters and forcings of the model run.

Return type:

Tuple[Float[Array, 'nz+1'], Float[Array, 'nz']]

Returns:

• bvf (float Array of shape (nz+1)) – Brunt–Väisälä frequency squared \(N^2\) on cell interfaces \(\left[\text s^{-2}\right]\).

• rho (Float[Array, ‘nz’]) – Density anomaly \(\rho\) on cell interfaces \(\left[\text {kg} \cdot \text m^{-3}\right]\)

Raises:

ValueError – If the value of case.eos_tracers is not one of {‘t’, ‘s’, ‘ts’, ‘b’}.

compute_shear(state, u_np1, v_np1)

Compute shear production term for TKE equation.

The prognostic equations are

\(S_h^2 = \partial_Z U^n \cdot \partial_z U^{n+1/2}\)

where \(U^{n+1/2}\) is the mean between \(U^n\) and \(U^{n+1}\).

Parameters:

• state (State) – Current state of the water column.

• u_np1 (float Array of shape (nz)) – Zonal velocity on the center of the cells at the next time step \(\left[\text m \cdot \text s^{-1}\right]\).

• v_np1 (float Array of shape (nz)) – Meridional velocity on the center of the cells at the next time step \(\left[\text m \cdot \text s^{-1}\right]\).

Returns:

shear2 – Shear production squared \(S_h^2\) on cell interfaces \(\left[\text m ^2 \cdot \text s ^{-3}\right]\).

Return type:

float Array of shape (nz+1)

compute_tke_eps_bc(tke, hz, keps_params, case_tracable)

Compute top and bottom boundary conditions for TKE and \(\varepsilon\).

Parameters:

• tke (float Array of shape (nz+1)) – Turbulent kinetic energy (TKE) denoted \(k\) on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-2}\right]\).

• hz (float Array of shape (nz)) – Thickness of cells from deepest to shallowest \([\text m]\).

• keps_params (KepsParameters) – Values of the parameters used by \(k-\varepsilon\).

• case_tracable (CaseTracable) – Physical parameters and forcings of the model run.

Return type:

Tuple[Float[Array, '1'], Float[Array, '1'], Float[Array, '1'], Float[Array, '1']]

Returns:

• tke_sfc_bc (float) – TKE value for surface boundary condition (Dirichlet \(\left[\text m ^2 \cdot \text s ^{-2}\right]\) or Neumann \(\left[\text m ^3 \cdot \text s ^{-3}\right]\)).

• tke_btm_bc (float) – TKE value for bottom boundary condition (Dirichlet \(\left[\text m ^2 \cdot \text s ^{-2}\right]\) or Neumann \(\left[\text m ^3 \cdot \text s ^{-3}\right]\)).

• eps_sfc_bc (float) – \(\varepsilon\) value for surface boundary condition (Dirichlet \(\left[\text m ^2 \cdot \text s ^{-3}\right]\) or Neumann \(\left[\text m ^3 \cdot \text s ^{-4}\right]\)).

• eps_btm_bc (float) – \(\varepsilon\) value for surface boundary condition (Dirichlet \(\left[\text m ^2 \cdot \text s ^{-3}\right]\) or Neumann \(\left[\text m ^3 \cdot \text s ^{-4}\right]\)).

Note

The kind of boundary conditions between Neumann and Dirichlet are register in the parameters keps_params.

advance_turb(akt, akv, tke, tke_np1, eps, c_mu, c_mu_prim, bvf, shear2, hz, dt, tke_sfc_bc, tke_btm_bc, eps_sfc_bc, eps_btm_bc, keps_params, do_tke)

Integrate TKE or \(\varepsilon\) quantities.

First the shear and buoyancy production are computed, then they are used in the building of the tridiagonal problem, the boundary conditions are then added and finally the tridiagonal problem is solved.

Parameters:

• akt (float Array of shape (nz+1)) – Current eddy-diffusivity on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-1}\right]\).

• akv (float Array of shape (nz+1)) – Current eddy-viscosity on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-1}\right]\).

• tke (float Array of shape (nz+1)) – Current turbulent kinetic energy (TKE) denoted \(k\) on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-2}\right]\).

• tke_np1 (float Array of shape (nz+1)) – Turbulent kinetic energy (TKE) denoted \(k\) on the interfaces of the cells at next step (usefull only for \(\varepsilon\) integration) \(\left[\text m ^2 \cdot \text s ^{-2}\right]\).

• eps (float Array of shape (nz+1)) – Current TKE dissipation denoted \(\varepsilon\) on the interfaces of the cells \(\left[\text m ^2 \cdot \text s ^{-3}\right]\).

• c_mu (float Array of shape (nz+1)) – Current \(c_\mu\) in GLS formalisim on the interfaces of the cells \([\text{dimensionless}]\).

• c_mu_prim (float Array of shape (nz+1)) – Current \(c_\mu'\) in GLS formalisim on the interfaces of the cells \([\text{dimensionless}]\).

• bvf (float Array of shape (nz+1)) – Current Brunt–Väisälä frequency squared \(N^2\) on cell interfaces \(\left[\text s^{-2}\right]\).

• shear2 (float Array of shape (nz+1)) – Current shear production squared \(S_h^2\) on cell interfaces \(\left[\text m ^2 \cdot \text s ^{-3}\right]\).

• hz (float Array of shape (nz)) – Thickness of cells from deepest to shallowest \([\text m]\).

• dt (float) – Time-step of the forward model \([\text s]\).

• tke_sfc_bc (float) – TKE value for surface boundary condition (Dirichlet \(\left[\text m ^2 \cdot \text s ^{-2}\right]\) or Neumann \(\left[\text m ^3 \cdot \text s ^{-3}\right]\)).

• tke_btm_bc (float) – TKE value for bottom boundary condition (Dirichlet \(\left[\text m ^2 \cdot \text s ^{-2}\right]\) or Neumann \(\left[\text m ^3 \cdot \text s ^{-3}\right]\)).

• eps_sfc_bc (float) – \(\varepsilon\) value for surface boundary condition (Dirichlet \(\left[\text m ^2 \cdot \text s ^{-3}\right]\) or Neumann \(\left[\text m ^3 \cdot \text s ^{-4}\right]\)).

• eps_btm_bc (float) – \(\varepsilon\) value for surface boundary condition (Dirichlet \(\left[\text m ^2 \cdot \text s ^{-3}\right]\) or Neumann \(\left[\text m ^3 \cdot \text s ^{-4}\right]\)).

• keps_params (KepsParameters) – Values of the parameters used by \(k-\varepsilon\).

• do_tke (bool) – If True solve the equation for TKE, else for \(\varepsilon\).

Returns:

vec – TKE or \(\varepsilon\) at next step (depending on do_tke).

Return type:

float Array of shape (nz+1)

compute_diag(tke, eps, bvf, shear2, keps_params)

Computes the diagnostic variables of \(k-\varepsilon\) closure.

This function first apply the Galperin limitation, then it computes \(c_\mu'\) and \(c_\mu\) with the stability function, and finally it computes the eddy-diffusivity and viscosity with these variables.

Parameters:

• tke (float Array of shape (nz+1)) – Turbulent kinetic energy (TKE) denoted \(k\) on the interfaces of the cells at next step \(\left[\text m ^2 \cdot \text s ^{-2}\right]\).

• eps (float Array of shape (nz+1)) – TKE dissipation denoted \(\varepsilon\) on the interfaces of the cells at next step \(\left[\text m ^2 \cdot \text s ^{-3}\right]\).

• bvf (float Array of shape (nz+1)) – Current Brunt–Väisälä frequency squared \(N^2\) on cell interfaces \(\left[\text s^{-2}\right]\).

• shear2 (float Array of shape (nz+1)) – Current shear production squared \(S_h^2\) on cell interfaces \(\left[\text m ^2 \cdot \text s ^{-3}\right]\).

• keps_params (KepsParameters) – Values of the parameters used by \(k-\varepsilon\).

Return type:

Tuple[Float[Array, 'nz+1'], Float[Array, 'nz+1'], Float[Array, 'nz+1'], Float[Array, 'nz+1'], Float[Array, 'nz+1']]

Returns:

• akt (float Array of shape (nz+1)) – Eddy-diffusivity on the interfaces of the cells at next step \(\left[\text m ^2 \cdot \text s ^{-1}\right]\).

• akv (float Array of shape (nz+1)) – Eddy-viscosity on the interfaces of the cells at next step \(\left[\text m ^2 \cdot \text s ^{-1}\right]\).

• eps (float Array of shape (nz+1)) – TKE dissipation denoted \(\varepsilon\) on the interfaces of the cells at next step \(\left[\text m ^2 \cdot \text s ^{-3}\right]\).

• c_mu (float Array of shape (nz+1)) – \(c_\mu\) in GLS formalisim on the interfaces of the cells at next step \([\text{dimensionless}]\).

• c_mu_prim (float Array of shape (nz+1)) – \(c_\mu'\) in GLS formalisim on the interfaces of the cells at next step \([\text{dimensionless}]\).