Implementation of an Eddy-Diffusivity Mass Flux (EDMF) scheme

[lemarief]

An Eddy-Diffusivity (ED) Mass-Flux (MF) scheme, as described in Perrot and Lemarié [2025] has been implemented in GOTM. The ED component is inherited from the 1-equation TKE scheme used in NEMO (Madec and al. [2024]) and inspired by Gaspar et al. [1990].

1. Eddy-Diffusivity component: NEMO 1 equation TKE scheme

For activating this option in the gotm.yaml file, the following options must be chosen:

turbulence :
   turb_method : first_order
   tke_method : tke
   len_scale_method : Gaspar_nemo
   stab_method : Prandtl_nemo

No other parameters are required because the internal parameter values of this scheme are fixed in the code to the standard NEMO values. The boundary conditions for the tke equation are chosen in the same way as for any other
turbulence scheme as the standard tkeeq subroutine is used. Most of the modifications were implemented in the file turbulence.F90 and concern the following elements:

integer , parameter , public :: Prandtl_nemo =4
integer , parameter , public :: Gaspar_nemo =7

In the do_lengthscale subroutine, there is a new case

case ( Gaspar_nemo )
   call lengthscale_gaspar (nlev ,z0s ,z0b ,h,NN)

The same in the stabilityfunctions subroutine,

case ( Prandtl_nemo )
   call cmue_nemo ( nlev )

Two new files ( lengthscale_gaspar.F90 and cmue_nemo.F90) have been added in the turbulence directory to implement the new option.

2. Mass-Flux component

The mass-flux component was much more cumbersome to implement, particularly in its
“energetically-consistent” version Perrot and Lemarié [2025]. Within this framework, the turbulent fluxes are expressed with a non-local component, in the form: $\overline{w' X'} = - K_X \partial_z \overline{X} + a_p w_p ( X_p - \overline{X} )$ with $X=u,v,\Theta,S$, $a_p$ the fractional area occupied by convective plumes and $w_p$ their vertical velocities. A set of nonlinear ordinary differential equations must be solved to compute $X_p, a_p$, and $w_p$ at each time-step. Moreover, these additional mass-flux terms affect the energy transfers and result in two extra terms in the TKE equation, as well as a modification of the turbulent transport of TKE: $P_{\rm MF} = - a_p w_p \partial_z \mathbf{u}\cdot(u_p-u) $, $G_{\rm MF} = a_p w_p ( b(X_p)-b(X) )$, and ${\cal D}_{\rm MF} = a_p w_p \left( k_p - k + \frac{1}{2} | \mathbf{u}_p - \overline{\mathbf{u}} |^2 \right)$.
Activating and configuring the 'mass flux’ option in the gotm.yaml file corresponds to the following options:

turbulence:
   MFconvec:
      compute_massflux: true
      massflux_on_dynamics: false
      massflux_energy: false
      mf_ap0: 0.2 
      mf_wp0: -0.5e-8 
      mf_Cent: 0.9 
      mf_Cdet: 1.8 
      mf_aa: 1.0 
      mf_bb: 1.0 
      mf_bp: 1.25 
      mf_uv: 0.25
      mf_dbkg: 1.25 

The compute_massflux logical controls the computation of the plume variables $a_p$ (fractional area), $w_p$ (plume vertical velocity), $T_p$ (plume temperature) and $S_p$ (plume salinity), and the application of the mass flux to the tracer equations. In
case massflux_on_dynamics is activated, $u_p$ and $v_p$ (plume horizontal velocities) are
computed, and a mass flux is applied to the u and v equations. Finally, when the logical
variable massflux_energy is set to true, additional TKE production terms $P_{\rm MF}$, $G_{\rm MF}$ are accounted for as well as ${\cal D}_{\rm MF}$.

3. Numerical validation for the implementation

The implementation has been validated using 2 idealized cases (free_convection and wind_convection) and 1 realistic case (asics_med) all presented in Perrot and Lemarié [2025]. The cases free_convection and wind_convection will be soon submitted as a pull request in the GOTM cases, and the realistic case is currently being integrated (see Bruno Deremble). Results obtained within GOTM are similar to those reported in Perrot and Lemarié [2025].

[bolding]

That is great with input from ' the outside' :-)

I think Lars is the better to judge the scientific contribution - but personally I'm not fond of a naming scheme using a specific model - nemo - is it not possible to make the attribution to a person?

Karsten

[lemarief]

We can refer to Madec rather than NEMO if needed; the appropriate reference would then be the NEMO documentation, as no publication specifically focused on this scheme (and how it differs from Gaspar et al. (1990) and Blanke & Delecluse (1993) ) has been published. Gurvan Madec has been the main contributor to the two specific components added in GOTM: the mixing lengths computation and limitation as well as the formulation of the turbulent Prandtl number.

Florian,

[lemarief]

More details about the mass-flux term implementation:

- Change in gotm.F90
After the vertical diffusion is applied and before computing the production term by shear,
the computation of mass flux variables (do_massflux) and the application to the meanflow temperature/salinity (and optionally horizontal momentum) equations is performed (uvTSequation_mf)

   if(compute_massflux) then
      call do_density(nlev,S,T,-z,-zi)
      call do_massflux( nlev,dt,u,v,T,S,rho,h )             ! compute mass flux quantities from the meanflow variables
      call uvTSequation_mf( nlev,dt,Fmass,u_p,v_p,T_p,S_p,Pmf,mf_nsub,massflux_on_dynamics, massflux_energy ) ! add the mass flux term to the meanflow equation 
   endif 

- Change in turbulence.F90

Additional arrays are declared

!  plumes properties for convective mass flux scheme 
   REALTYPE, public, dimension(:), allocatable   :: a_p,w_p,Fmass
   REALTYPE, public, dimension(:), allocatable   :: T_p,S_p,wk_mf 
   REALTYPE, public, dimension(:), allocatable   :: normVp,tke_p
   REALTYPE, public, dimension(:), allocatable   :: u_p,v_p,EmD

In the do_turbulence subroutine a mass-flux term on tke can be optionally applied

   case (first_order)
      call production(nlev,NN,SS, xP=xP, SSCSTK=SSCSTK, SSSTK=SSSTK)
      call alpha_mnb(nlev,NN,SS)
      call stabilityfunctions(nlev)
      call do_tke(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
      if(compute_massflux .and. massflux_energy) then 
         call tkeeq_mf(nlev,dt,h)
      endif
      call do_lengthscale(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS)
      call kolpran(nlev)
      call internal_wave(nlev,NN,SS)

a new do_massflux subroutine has been added, this subroutine is responsible for the computation of plume properties from meanflow variables

   call massfluxeq(nlev,dt,umean,vmean,tmean,smean,Rmean,h)
   if( massflux_energy ) call tke_transport_mf(nlev,umean,vmean,h)

- Change in tkeeq.F90
Finally, the TKE equation is adapted when energy consistency (and thus additional TKE production terms) is activated.

   do i=1,nlev-1

!     compute diffusivity
      avh(i) = num(i)/sig_k

!     compute production terms in k-equation
      prod     = P(i) + Px(i) + PSTK(i)
      buoyan   = B(i) + Bmf(i) + Pmf(i) + wk_mf(i)
      diss     = eps(i)


      if (prod+buoyan.gt.0) then
         Qsour(i)  = prod+buoyan
         Lsour(i) = -diss/tke(i)
      else
         Qsour(i)  = prod
         Lsour(i) = -(diss-buoyan)/tke(i)
      end if

   end do

the extra terms Bmf, Pmf, and wk_mf are added to the buoyan component as they are not necessarily positive-definite.