Flake methane scheme is supplemented with treatment of the case when the mixed...

Flake methane scheme is supplemented with treatment of the case when the mixed layer touches the lake bottom

source/Flake/flake_methane.f90

@@ -6,6 +6,7 @@ implicit none
 
 real(kind=ireals), parameter :: cw    = 3990.          !> Specific heat of water,  J/(kg*K)
 real(kind=ireals), parameter :: row0 = 1.E+3           !> Reference water density, kg/m**3
+real(kind=ireals), parameter :: niu_wat = 1.307d-6     !> Molecular viscosity of water at T = 10 C, m**2/s
 real(kind=ireals), parameter :: g = 9.80665            !> Reference acceleration due to gravity, m/s**2
 real(kind=ireals), parameter :: Kelvin0 = 273.16       !> Temperature of melting point, Kelvin
 real(kind=ireals), parameter :: pres_ref = 1013.5*100. !> Atmospheric reference pressure, Pa
@@ -28,14 +29,14 @@ real(kind=ireals), parameter :: Henry_temp_dep_o2 = 1.7d+3  !> The temperature d
 
 real(kind=ireals), parameter :: molm3tomgl_o2 = 32.
 real(kind=ireals), parameter :: molm3tomgl_ch4 = 16.
-real(kind=ireals), parameter :: k_ch4 = 0.6 / molm3tomgl_ch4 !(Liikanen et al. 2002; Lofton et al. 2013)
-                                                !0.44 ! Michaelis constant for methane in methane oxidation,
-                                                ! after Arah&Stephen (1998), Nedwell&Watson (1995) 
-real(kind=ireals), parameter :: k_o2 = 0.672 / molm3tomgl_o2 ! (Lidstrom and Somers, 1984)
-                                                !0.33  ! Michaelis constant for oxygen in methane oxidation, 
-                                                ! after Arah&Stephen (1998), Nedwell&Watson (1995)
-real(kind=ireals), parameter :: Vmaxw = 1.d-1/86400. ! reaction potential in oxygen-saturated Michaelis-Menten kinetics, 
-                                        ! after (Liikanen et al., 2002)
+real(kind=ireals), parameter :: k_ch4 = 0.6 / molm3tomgl_ch4 !>(Liikanen et al. 2002; Lofton et al. 2013)
+                                                             !!0.44 ! Michaelis constant for methane in methane oxidation,
+                                                             !! after Arah&Stephen (1998), Nedwell&Watson (1995) 
+real(kind=ireals), parameter :: k_o2 = 0.672 / molm3tomgl_o2 !> (Lidstrom and Somers, 1984)
+                                                             !!0.33  ! Michaelis constant for oxygen in methane oxidation, 
+                                                             !! after Arah&Stephen (1998), Nedwell&Watson (1995)
+real(kind=ireals), parameter :: Vmaxw = 1.d-1/86400. !> Reaction potential in oxygen-saturated Michaelis-Menten kinetics, 
+                                                     !! after (Liikanen et al., 2002)
 
 
 contains
@@ -47,12 +48,13 @@ implicit none
 real(kind=ireals) :: HENRY_CONST
 
 ! Input variables
-real(kind=ireals), intent(in) :: henry_const0 ! Henry constant at the reference temperature, mol/(m**3*Pa)
-real(kind=ireals), intent(in) :: temp_dep ! Temperature dependence (enthalpy solution devided by universal gas constant), K
-real(kind=ireals), intent(in) :: temp_ref ! Reference temperature, K
-real(kind=ireals), intent(in) :: temp ! Temperature, K
+real(kind=ireals), intent(in) :: henry_const0 !> Henry constant at the reference temperature, mol/(m**3*Pa)
+real(kind=ireals), intent(in) :: temp_dep     !> Temperature dependence (enthalpy solution devided by universal gas constant), K
+real(kind=ireals), intent(in) :: temp_ref     !> Reference temperature, K
+real(kind=ireals), intent(in) :: temp         !> Temperature, K
 ! real(kind=ireals), intent(in) :: radius ! Curvature radius, m, positive for drops, negative for bubbles and zero for flat surface
 
+
 HENRY_CONST = henry_const0*exp(temp_dep*(1./temp-1./temp_ref))
 ! Effect of bubble surface curvatue on saturation pressure
 !  if (radius /= 0.d0) HENRY_CONST = HENRY_CONST* &
@@ -101,7 +103,26 @@ k = sqrt((C1*U)**2 + (C2*wstar)**2)*SCHMIDT_NUMBER_METHANE(temp-Kelvin0)**(-0.5)
 END FUNCTION PISTON_VELOCITY_HEISKANEN
 
 
-SUBROUTINE METHANE_MAIN(H,h_ML,h_s,wstar,U,pres,T_ML,Tb,C_ML,F_difsurf,F_bubble)
+!> Function calculates the thickness of the bottom viscous sublayer
+FUNCTION THICKNESS_VISC_SUBLAYER(ustar,wstar) result(hvisc)
+
+implicit none
+
+!Input variables
+real(kind=ireals), intent(in) :: wstar !> Convective Deardorff velocity scale in the mixed layer, m/s
+real(kind=ireals), intent(in) :: ustar !> Friction velocity scale in the mixed layer, m/s
+
+!Local variables
+real(kind=ireals), parameter :: C  = 8.5 !> 5 -- 11.6 (Wengrove and Foster, GRL 2014, and references therein)
+real(kind=ireals), parameter :: C1 = 1.  !> @todo: to be checked in literature
+real(kind=ireals) :: hvisc
+
+hvisc = C*niu_wat/sqrt(ustar*ustar + C1 * wstar*wstar)
+
+END FUNCTION THICKNESS_VISC_SUBLAYER
+
+
+SUBROUTINE METHANE_MAIN(H,h_ML,h_s,wstar,ustar,U,pres,T_ML,Tb,C_ML,F_difsurf,F_bubble)
 
 implicit none
 
@@ -110,6 +131,7 @@ real(kind=ireals), intent(in) :: H     !> Lake depth, m
 real(kind=ireals), intent(in) :: h_ML  !> Mixed-layer depth, m
 real(kind=ireals), intent(in) :: h_s   !> The thickness of bottom sediments layer, m
 real(kind=ireals), intent(in) :: wstar !> Convective Deardorff velocity scale in the mixed layer, m/s
+real(kind=ireals), intent(in) :: ustar !> Friction velocity scale in the mixed layer, m/s
 real(kind=ireals), intent(in) :: U     !> Wind speed, m/s
 real(kind=ireals), intent(in) :: pres  !> Atmospheric pressure, Pa
 real(kind=ireals), intent(in) :: T_ML  !> Mixed-layer temperature, Kelvin
@@ -120,7 +142,7 @@ real(kind=ireals), intent(out) :: C_ML      !> Methane concentration in the mixe
 real(kind=ireals), intent(out) :: F_difsurf !> Surface diffusive methane flux to the atmosphere, mol/(m**3*s)
 real(kind=ireals), intent(out) :: F_bubble  !> Bubble methane flux to the atmosphere, mol/(m**3*s)
 
-real(kind=ireals) :: k, k_, sin1, cos1
+real(kind=ireals) :: k, k_, sin1, cos1, h_ML_eff, dh
 real(kind=ireals) :: gas_exch_const        !> Piston velocity for methane at the water surface, m/s
 real(kind=ireals) :: C_b                   !> Methane concentration at the lake bottom, mol/m**3
 real(kind=ireals) :: C_O2_ML               !> Oxygen concentration in the mixed layer, mol/m**3
@@ -135,6 +157,10 @@ real(kind=ireals) :: diff_hypo             !> Diffusion coefficient in hypolimni
 real(kind=ireals) :: meth_oxid_const_hypo  !> Methane oxidation rate constant in hypolimnion
 real(kind=ireals) :: meth_oxid_const_ML    !> Methane oxidation rate constant in the mixed layer
 
+! Limiting hypolimnion thickness from below by the bottom viscous sublayer depth
+dh = max (H - h_ML, THICKNESS_VISC_SUBLAYER(ustar,wstar))
+h_ML_eff = H - dh
+
 ! Calculation of methane concentration in the mixed layer, equilibrated with atmospheric concentration
 CH4ppb = 1.8E+3
 CH4pres = CH4ppb*1.E-9*pres
@@ -165,11 +191,11 @@ meth_oxid_const_hypo = Vmaxw/k_ch4*C_O2_hypo/(k_o2 + C_O2_hypo)
 
 k = sqrt(meth_oxid_const_hypo/diff_hypo)
 k_ = diff_hypo*k
-sin1 = sinh(k*(H - h_ML))
-cos1 = cosh(k*(H - h_ML))
+sin1 = sinh(k*(H - h_ML_eff))
+cos1 = cosh(k*(H - h_ML_eff))
 
 C_ML = (k_ * C_b + gas_exch_const * meth_atmo * sin1) / &
-     & ( (gas_exch_const + meth_oxid_const_ML * h_ML * sin1 ) + k_* cos1 )
+     & ( (gas_exch_const + meth_oxid_const_ML * h_ML_eff * sin1 ) + k_* cos1 )
 
 Prod = Prod0 * q10 ** (0.1*(Tb - Kelvin0)) / alpha0 * &
        & (1. - exp(-alpha0*h_s))

source/Flake/src_flake_interface_1D.F90

@@ -393,7 +393,7 @@ CALL flake_driver ( depth_w, depth_bs, T_bs, par_Coriolis,         &
 !  Update methane concentration and fluxes (Victor Stepanenko, 2017/07/31/)
 !------------------------------------------------------------------------------
 
-call METHANE_MAIN(depth_w, h_ML_n_flk, depth_bs, w_star_sfc_flk, &
+call METHANE_MAIN(depth_w, h_ML_n_flk, depth_bs, w_star_sfc_flk, u_star_w_flk, &
                   U_a_in, P_a_in, T_wML_n_flk, T_bot_n_flk, C_ML, F_difsurf, F_bubble)