\[\newcommand{\paramT}[1]{ \text{#1}}
\newcommand{\hyperparamT}[1]{\text{#1}}
\newcommand{\simparamT}[1]{ \text{#1}}
\newcommand{\exp}[1]{\mathrm{exp}\left(#1\right)}
\newcommand{\atan}[1]{\mathrm{atan}\left(#1\right)}
\newcommand{\sign}[1]{\mathrm{sign}\left(#1\right)}
\newcommand{\erf}[1]{\mathrm{erf}\left(#1\right)}
\newcommand{\erfinv}[1]{\mathrm{erfinv}\left(#1\right)}
\newcommand{\param}[1]{ #1}
\newcommand{\hyperparam}[1]{#1}
\newcommand{\simparam}[1]{ #1}
\newcommand{\CROSS}{\times}
\newcommand{\GRAD}{\nabla}
\newcommand{\DOT}{\bullet}
\newcommand{\PD}{\partial}
\newcommand{\PDFz}{\frac{\PD}{\PD z}}
\newcommand{\DM}[1]{\langle #1 \rangle}
\newcommand{\iEnv}{e}
\newcommand{\SD}[2]{{\overline{#1}}_{#2}}
\newcommand{\SDi}[1]{{\SD{#1}{i}}}
\newcommand{\SDj}[1]{{\SD{#1}{j}}}
\newcommand{\SDe}[1]{{\SD{#1}{\iEnv}}}
\newcommand{\SDiog}[2]{#1_{#2}}
\newcommand{\SDio}[1]{{\SDiog{#1}{i}}}
\newcommand{\SDjo}[1]{{\SDiog{#1}{j}}}
\newcommand{\SDeo}[1]{{\SDiog{#1}{\iEnv}}}
\newcommand{\aSD}[2]{{#1}_{#2}}
\newcommand{\aSDi}[1]{\aSD{#1}{i}}
\newcommand{\aSDj}[1]{\aSD{#1}{j}}
\newcommand{\aSDe}[1]{\aSD{#1}{\iEnv}}
\newcommand{\otherDefs}{where additional variable definitions are in:}
\newcommand{\IntraCVSDi}[2]{\overline{{#1}_{i }'{#2}_{i }'}}
\newcommand{\IntraCVSDj}[2]{\overline{{#1}_{j }'{#2}_{j }'}}
\newcommand{\IntraCVSDe}[2]{\overline{{#1}_{\iEnv{}}'{#2}_{\iEnv{}}'}}
\newcommand{\InterCVSDi}[2]{\overline{{#1}_{i }'}~\overline{{#2}_{i }'}}
\newcommand{\InterCVSDj}[2]{\overline{{#1}_{j }'}~\overline{{#2}_{j }'}}
\newcommand{\InterCVSDe}[2]{\overline{{#1}_{\iEnv{}}'}~\overline{{#2}_{\iEnv{}}'}}
\newcommand{\TCV}[2]{\langle {#1}^*{#2}^* \rangle}
\newcommand{\BC}[1]{{#1|_{z_{min}}}}
\newcommand{\BCT}[1]{{#1|_{z_{max}}}}
\newcommand{\BCB}[1]{{#1|_{z_{min}}}}
\newcommand{\BCG}[1]{{#1|_{z_{boundary}}}}
\newcommand{\Km}{K^m}
\newcommand{\Kh}{K^h}
\newcommand{\TEquilib}{T_{\mathrm{iterated}}}
\newcommand{\PhasePartition}{q}
\newcommand{\ExnerD}{\Pi_{dry}}
\newcommand{\ExnerM}{\Pi_{moist}}
\newcommand{\WindSpeed}{|u|}
\newcommand{\LayerThickness}{\param{\Delta z}}
\newcommand{\SurfaceRoughness}[1]{\param{z_{0#1}}}
\newcommand{\SensibleSurfaceHeatFlux}{F_{\mathrm{sensible}''}}
\newcommand{\LatentSurfaceHeatFlux}{F_{\mathrm{latent}''}}
\newcommand{\FrictionVelocity}{u_*}
\newcommand{\Buoyancy}{b}
\newcommand{\BuoyancyGrad}{\PD_z \Buoyancy}
\newcommand{\BuoyancyFlux}{\IntraCVSDi{w}{b}}
\newcommand{\TemperatureScale}{\theta_*}
\newcommand{\SurfaceMomentumFlux}{\BC{\overline{w'u'}}}
\newcommand{\SurfaceHeatFlux}{\BC{\overline{w'\theta'}}}
\newcommand{\SurfaceBuoyancyFlux}{\BC{\IntraCVSDi{w}{\theta}}}
\newcommand{\ConvectiveVelocity}{{w_*}} % Convective velocity near the surface
\newcommand{\InversionHeight}{{z_*}}
\newcommand{\MOLen}{\Lambda_{M-O}}
\newcommand{\zLL}{\param{z_{||}}} % z at the first surface level (we should make this grid-independent)
\newcommand{\qt}{q_{\mathrm{tot}}}
\newcommand{\qr}{q_{\mathrm{rain}}}
\newcommand{\ql}{q_{\mathrm{liq}}}
\newcommand{\qi}{q_{\mathrm{ice}}}
\newcommand{\qv}{q_{\mathrm{vap}}}
\newcommand{\qvsat}{q_{\mathrm{vap}}^*}
\newcommand{\pvsat}{p_{\mathrm{vap}}^*}
\newcommand{\qc}{q_{\mathrm{con}}}
\newcommand{\ThetaVap}{{\theta_{\mathrm{vap}}}}
\newcommand{\ThetaVirt}{{\theta_{\mathrm{virt}}}}
\newcommand{\ThetaRho}{{\theta_{\rho}}}
\newcommand{\ThetaLiq}{{\theta_{\mathrm{liq}}}}
\newcommand{\ThetaLiqIce}{{\theta_{\mathrm{liq-ice}}}}
\newcommand{\ThetaLiqIceSat}{{\theta^*_{\mathrm{liq-ice}}}}
\newcommand{\ThetaDry}{{\theta_{\mathrm{dry}}}}
\newcommand{\TDry}{{T_{dry}}}
\newcommand{\eint}{e_{\mathrm{int}}}
\newcommand{\etot}{e_{\mathrm{tot}}}
\newcommand{\TRef}{{T}_0}
\newcommand{\alphaRef}{{\alpha}_0}
\newcommand{\rhoRef}{{\rho}_0}
\newcommand{\pRef}{{p}_0}
\newcommand{\Heaviside}{\mathcal H}
\newcommand{\alphaLL}{\alphaRef|_{\zLL}}
\newcommand{\uH}{\simparam{\mathbf{u}_h}}
\newcommand{\CoriolisParam}{\hyperparam{\mathrm{coriolis\_param}}}
\newcommand{\SubsidenceParam}{\hyperparam{\mathrm{subsidence}}}
\newcommand{\betaM}{\hyperparam{\beta_m}}
\newcommand{\betaH}{\hyperparam{\beta_h}}
\newcommand{\gammaM}{\hyperparam{\gamma_m}}
\newcommand{\gammaH}{\hyperparam{\gamma_h}}
\newcommand{\PTilde}{\param{\tilde{p}}}
\newcommand{\VKConst}{\param{\kappa_{\mathrm{Von-Karman}}}}
\newcommand{\Nsd}{\hyperparam{N_{sd}}}
\newcommand{\grav}{\param{g}}
\newcommand{\TZero}{\param{T_{0}}}
\newcommand{\RefHintV}{\param{{\eint}_{v,0}}}
\newcommand{\RefHintI}{\param{{\eint}_{i,0}}}
\newcommand{\EpsDV}{\param{\varepsilon_{dv}}}
\newcommand{\EpsVD}{\param{\varepsilon_{vd}}}
\newcommand{\Rm}{R_m}
\newcommand{\Cpm}{c_{pm}}
\newcommand{\Cvm}{c_{vm}}
\newcommand{\Rd}{\param{R_d}}
\newcommand{\Rv}{\param{R_v}}
\newcommand{\Cp}[1]{\param{c_{p#1}}}
\newcommand{\Cv}[1]{\param{c_{v#1}}}
\newcommand{\Cvd}{\Cv{d}}
\newcommand{\Cvv}{\Cv{v}}
\newcommand{\Cvl}{\Cv{l}}
\newcommand{\Cvi}{\Cv{i}}
\newcommand{\DeltaCp}{\param{\Delta c_p}}
\newcommand{\TTriple}{\param{T_{\mathrm{tr}}}}
\newcommand{\PTriple}{\param{p_{\mathrm{tr}}}}
\newcommand{\TFreeze}{\param{T_{\mathrm{freeze}}}}
\newcommand{\RefLHv}{\param{L_{v,0}}}
\newcommand{\RefLHs}{\param{L_{s,0}}}
\newcommand{\RefLHf}{\param{L_{f,0}}}
\newcommand{\LatentHeatV}[1]{L_{vap}(#1)}
\newcommand{\LatentHeatS}[1]{L_{sub}(#1)}
\newcommand{\LatentHeatF}[1]{L_{fus}(#1)}\]
This document describes the Eddy-Diffusivity Mass-Flux (EDMF) parameterization for subgrid-scale (SGS) turbulence and convection. The model equations and rationale are based on: Tan et al. (2018); Cohen et al. (2020); Lopez-Gomez et al. (2020). The key predictands of the EDMF parameterization are the SGS vertical fluxes of heat, moisture and momentum, and the cloud fraction in the host model grid. The EDMF parameterization divides the host model's grid into $N \ge 2$ subdomains that represent an isotropic turbulent enviroment and coherent updraft(s) and or downdraft(s). The parameterization solves prognostic equations for first and second moments in the subdomains to provide both the abovementioned SGS vertical fluxes and cloud fraction in the host model grid.
In this document, color-coding is used to indicate:
\[\paramT{Constant parameters that are fixed in space and time (e.g., those defined in CLIMAParameters.jl)}\]
\[\simparamT{Single column (SC) inputs (e.g., variables that are fed into the SC model from the dynamical core (e.g., horizontal velocity))}\]
\[\hyperparamT{Tunable hyper-parameters that will need to be changeable, but will only include single numbers (e.g., Float64)}\]
The EDMF is 1D model in $z$ that solves for the statistical distribution in the grid box of the host model. As such in the EDMF there is no spatial discretization in the horizontal directions ($x$ and $y$), the horizontal space is broken into $\Nsd$ ($\sim$ 1-5) "subdomains" (SDs), denoted by subscript $i$, where $1 \le i \le \Nsd$. One of the subdomains, the "environment", is treated different compared to others, termed "updrafts" (and or "downdrafts"). This environment subdomain is denoted with a special index $\iEnv{}$ (which we usually set to 0). For dummy variables $\phi$ and $\psi$, we use several domain and SD representations of interest:
\[\begin{align}
\SDi{\phi} \quad & \text{horizontal mean of variable $\phi$ over subdomain $i$}, \\
\SDi{\phi}' = \phi_i - \SDi{\phi} \quad & \text{fluctuations of $\phi$ about the subdomain mean}, \\
\IntraCVSDi{\phi}{\psi} \quad & \text{intra subdomain covariance}, \\
\DM{\phi} = \sum_i \aSDi{a} \SDi{\phi} \quad & \text{horizontal mean of $\phi$ over the total domain}, \\
\SDi{\phi}^* = \SDi{\phi} - \DM{\phi} \quad & \text{difference between subdomain & domain means}, \\
\InterCVSDi{\phi}{\psi} \quad & \text{inter subdomain covariance among subdomain means}, \\
\phi^* = \phi - \DM{\phi} \quad & \text{difference between subdomain & domain means}, \\
\TCV{\phi}{\psi} = \sum_{\forall i} a_i \IntraCVSDi{\phi}{\psi} +
\sum_{\forall i} \sum_{\forall j} \aSDi{a} \aSDj{a} \SDi{\phi}(\SDi{\psi} - \SDj{\psi}) \quad & \text{total covariance}.
\end{align}\]
Here, $\SDi{\phi}$ and $\SDi{\psi}$ are a dummy variables for the following 7 unknowns:
\[\begin{align}
\SDi{w} & \quad \text{vertical velocity}, \\
\SDi{\eint} & \quad \text{internal energy}, \\
\SDi{\qt} & \quad \text{total water specific humidity}, \\
\SDi{TKE} & \quad \text{turbulent kinetic energy ($0.5(\IntraCVSDi{u}{u}+\IntraCVSDi{v}{v}+\IntraCVSDi{w}{w})$)}, \\
\IntraCVSDi{\eint}{\eint} & \quad \text{intra subdomain covariance of $\eint'$ and $\eint'$}, \\
\IntraCVSDi{\qt}{\qt} & \quad \text{intra subdomain covariance of $\qt'$ and $\qt'$}, \\
\IntraCVSDi{\eint}{\qt} & \quad \text{intra subdomain covariance of $\eint'$ and $\qt'$}.
\end{align}\]
From the large-scale model perspective, $\DM{\phi}$ represents the resolved grid mean, and $\TCV{\phi}{\psi}$ represents the SGS fluxes and (co)-variances of scalars that need to be parameterized. Equations in the following sections, \eqref{eq:AreaFracGov}, \eqref{eq:1stMoment} and \eqref{eq:2ndMoment}, are solved on $z_{min} \le z \le z_{max}$ and $t \ge 0$. There are $8 \Nsd$ equations in total.
The EDMF model can be used in the context of a stand-alone single column, or integrated with a dynamical core. Either way, the EDMF model relies on grid mean variables, which may be prescribed or solved for. Taking an area fraction-weighted average of the subdomain equations yields the domain-averaged equations (which should be consistent with variables in the dynamical core).
The domain-averaged equations for $\DM{\phi} \in [w, \qt, \eint, \uH]$ are:
\[\begin{align}
\PD_t (\rhoRef{} \DM{\phi})
+ \PD_z (\rhoRef{} \DM{w} \DM{\phi})
+ \nabla_h \DOT \left( \rhoRef{} \DM{\phi} \otimes \DM{\phi} \right)
= \\
\DM{S}_{\text{diff}}^{\DM{\phi}}
+ \DM{S}_{\text{press}}
+ \DM{S}_{\text{coriolis}}
+ \DM{S}_{\text{subsidence}},
\end{align}\]
where
\[\begin{align}
\DM{S}_{\text{diff}}^{\DM{\phi}} & = \PD_z (\rhoRef{} \aSDe{a} \SDe{\Km} \PD_z \DM{\phi}), \label{eq:gm_diffusion} \\
\DM{S}_{\text{diff}}^{w} & = \PD_z (\rhoRef{} \aSDe{a} \SDe{\Km} \PD_z \DM{w}) , \label{eq:gm_diffusion_w} \\
\DM{S}_{\text{press}} & = - \GRAD_h \DM{p}, \label{eq:gm_pressure} \\
\DM{S}_{\text{coriolis}} & = \CoriolisParam \DM{\phi} \CROSS \mathbf{k}, \label{eq:gm_coriolis} \\
\DM{S}_{\text{subsidence}} & = - \SubsidenceParam \GRAD \phi, \label{eq:gm_subsidence} \\
\end{align}\]
The EDMF equations take the form of advection-diffusion equations. The subdomain area fraction is given by the mass continuity equation in the $i$th subdomain ($\aSDi{a}$):
\[\begin{gather}
\PD_t (\rhoRef{} \aSDi{a})
+ \PD_z (\rhoRef{} \aSDi{a} \SDi{w})
+ \GRAD_h \DOT
(\rhoRef{} \aSDi{a} \DM{\uH})
=
\SDi{S}^a
, \quad i \ne \iEnv, \label{eq:AreaFracGov} \\
\aSDi{a} = 1 - \sum_{j\ne\iEnv} \aSDj{a}, \quad i = \iEnv, \label{eq:AreaFracConserve} \\
\qquad 0 < \aSDi{a} < 1. \label{eq:AreaFracConstraint}
\end{gather}\]
Here, $\rhoRef{}, \SDi{w}, \uH$ is fluid density, mean vertical velocity along $z$, and domain-mean of the horizontal velocity respectively. The area fraction constraints are necessary to ensure the system of equations is well-posed. All source terms ($\SDi{S}^a$) will be discussed in later sections.
The greater than zero constraint must be satisfied at every step of the solution process, since it is necessary to avoid division by zero in the mean field equations.
\[\begin{align}
\SDi{S}^a = \SDi{S_{\epsilon\delta}}^a.
\end{align}\]
We note that the net exchange is zero $\sum_i \SDi{S_{\epsilon\delta}}^a = 0$. Therefore, we may define the environment source term as the negative sum of all updraft source terms. The entrainment-detrainment source is:
\[\begin{align}
\SDi{S_{\epsilon\delta}}^a =
\begin{cases}
\rho a_i \SDi{w} \left( -\delta_i + \epsilon_{i} \right) & i \ne \iEnv{} \\
0 - \sum_{j \ne \iEnv{}} \SDj{S_{\epsilon\delta}}^a & i = \iEnv{} \\
\end{cases},
\end{align}\]
where additional variable definitions are in:
The 1st moment sub-domain equations are:
\[\begin{align}\label{eq:1stMoment}
\PD_t (\rhoRef{} \aSDi{a} \SDi{\phi})
+ \PD_z (\rhoRef{} \aSDi{a} \SDi{w} \SDi{\phi})
+ \GRAD_h \DOT
(\rhoRef{} \aSDi{a} \DM{\uH} \SDi{\phi})
=
\SDi{S}^\phi
, \quad \forall i. \\
\end{align}\]
Here, $\SDi{S}^{\phi}$ are source terms, including diffusion, and many other sub-grid-scale (SGS) physics. In general, $\SDi{S}^{\phi}$ and $\SDi{S}^{a}$ may depend on $\SDj{\phi}$ and or $\aSDj{a}$ for any $j$.
The source terms common to all unknowns are:
\[\begin{align}
\SDi{S}^{\phi} =
\SDi{S_{\epsilon\delta}}^{\phi}
+ \SDi{S_{\text{turb-transp}}}^{\phi}, \quad \forall \phi
\end{align}\]
Additional source terms exist in other equations:
\[\begin{align}
\SDi{S}^{w} &=
\SDi{S_{\epsilon\delta}}^w
+ \SDi{S_{\text{turb-transp}}}^w
+ \SDi{S_{\text{buoy}}}
+ \SDi{S_{\text{nh-press}}}
+ \SDi{S_{\text{coriolis}}}, \\
\SDi{S}^{\eint} &=
\SDi{S_{\epsilon\delta}}^{\eint}
+ \SDi{S_{\text{turb-transp}}}^{\eint}
+ \SDi{S_{\text{MP-MSS}}}^{\eint}
+ \SDi{S_{\text{rad}}}, \\
\SDi{S}^{\qt} &=
\SDi{S_{\epsilon\delta}}^{\qt}
+ \SDi{S_{\text{turb-transp}}}^{\qt}
+ \SDi{S_{\text{MP-MSS}}}^{\qt}.
\end{align}\]
Note: The sum of the total pressure and gravity are recast into the sum of the non-hydrostatic pressure and buoyancy sources.
\[\begin{align}
\SDi{S_{\epsilon\delta}}^{\phi} &=
\begin{cases}
\rhoRef{} a_i \SDi{w} \left( -\delta_i \SDi{\phi} + \epsilon_{i} \SDe{\phi} \right) & i \ne \iEnv{} \\
0 - \sum_{j \ne \iEnv{}} \SDj{S_{\epsilon\delta}}^{\phi} & i=\iEnv{} \\
\end{cases} \\
\SDi{S_{\text{turb-transp}}}^{\phi} & = -\PD_z (\rhoRef{} a_i \IntraCVSDi{w}{\phi}) \\
& = \PD_z (\rhoRef{} a_i \SDi{\Km} \PD_z \SDi{\phi}) \\
\SDi{S_{\text{nh-press}}} &= -\rhoRef{} \aSDi{a} \left( \alpha_b \SDi{b} + \alpha_d \frac{(\SDi{w} - \SDe{w}) | \SDi{w} - \SDe{w} | }{r_d \aSDi{a}^{1/2}} \right) \\
\alpha_b &= 1/3, \quad \alpha_d = 0.375, \quad r_d = 500 [m] \\
\SDi{S_{\text{buoy}}} &= \rhoRef{} \aSDi{a} \SDi{b} \\
\SDi{S_{\text{coriolis}}} & = f(\SDi{\mathbf{u}} - {\SDi{\mathbf{u}_{\text{geo-wind}}}}) \\
\SDi{S_{\text{rad}}} &= \left( \PD_t {\SDi{\eint}} \right)_{radiation} \\
\SDi{S_{\text{MP-MSS}}}^{\qt} & = \\
\SDi{S_{\text{MP-MSS}}}^{\eint} & = \\
\end{align}\]
where additional variable definitions are in:
The 2nd moment sub-domain equations are of the exact same form as the 1st moment equations (equation \eqref{eq:1stMoment}):
\[\begin{align}\label{eq:2ndMoment}
\PD_t (\rhoRef{} \aSDi{a} \SDi{\phi})
+ \PD_z (\rhoRef{} \aSDi{a} \SDi{w} \SDi{\phi})
+ \GRAD_h \DOT
(\rhoRef{} \aSDi{a} \DM{\uH} \SDi{\phi})
=
\SDi{S}^\phi
, \quad \forall i. \\
\end{align}\]
Here, $\SDi{S}^{\phi}$ are source terms, including diffusion, and many other sub-grid-scale (SGS) physics. In general, $\SDi{S}^{\phi}$ and $\SDi{S}^{a}$ may depend on $\SDj{\phi}$ and or $\aSDj{a}$ for any $j$.
The source terms common to all unknowns are:
\[\begin{align}
\SDi{S}^{\phi\psi} =
\SDi{S_{\epsilon\delta}}^{\phi\psi}
+ \SDi{S_{\text{x-grad flux}}}^{\phi\psi}
+ \SDi{S_{\text{turb-transp}}}^{\phi\psi}
\quad \forall \phi \psi,
\end{align}\]
Additional source terms exist in other equations:
\[\begin{align}
\SDi{S}^{TKE} &=
\SDi{S_{\epsilon\delta}}^{TKE}
+ \SDi{S_{\text{x-grad flux}}}^{TKE}
+ \SDi{S_{\text{turb-transp}}}^{TKE}
+ \SDi{S_{\text{dissip}}}
+ \SDi{S_{\text{press}}},
+ \SDi{S_{\text{buoyancy}}}, \\
\SDi{S}^{\phi\psi} &=
\SDi{S_{\epsilon\delta}}^{\phi\psi}
+ \SDi{S_{\text{x-grad flux}}}^{\phi\psi}
+ \SDi{S_{\text{turb-transp}}}^{\phi\psi}
+ \SDi{S_{\text{dissip}}}^{\phi\psi}
+ \SDi{S_{\text{MP-MSS}}}^{\phi\psi}.
\quad \phi\psi \in [\qt\qt, \eint\eint, \eint \qt].
\end{align}\]
\[\begin{align}
\SDi{S_{\epsilon\delta}}^{\phi\psi} &=
\begin{cases}
\rhoRef{} a_i \SDi{w} \left[ -\delta_i \IntraCVSDi{\phi}{\psi} + \epsilon_{i}
\left(
\IntraCVSDe{\phi}{\psi} + (\SDe{\phi} - \SDi{\phi})(\SDe{\psi} - \SDi{\psi})
\right) \right] & i \ne \iEnv \\
0 - \sum_{j\ne \iEnv} \SDj{S_{\epsilon\delta}}^{\phi\psi} & i=\iEnv \\
\end{cases} \\
\SDi{S_{\epsilon\delta}}^{TKE} &=
\begin{cases}
\rhoRef{} a_i \SDi{w} \left[ -\delta_i \SDi{TKE} + \epsilon_{i}
\left(
\SDe{TKE} + \frac{1}{2} (\SDe{w} - \SDi{w})^2
\right) \right] & i \ne \iEnv \\
0 - \sum_{j\ne \iEnv} \SDj{S_{\epsilon\delta}}^{TKE} & i=\iEnv \\
\end{cases} \\
\SDi{S_{\text{x-grad flux}}}^{\phi\psi}
& =
- \rhoRef{} a_i \IntraCVSDi{w}{\psi} \PD_z \SDi{\phi}
- \rhoRef{} a_i \IntraCVSDi{w}{\phi} \PD_z \SDi{\psi} \\
& =
2 \rhoRef{} a_i \SDi{\Kh} \PD_z \SDi{\psi} \PD_z \SDi{\phi} \\
\SDi{S_{\text{x-grad flux}}}^{TKE}
& =
\rhoRef{} a_i \SDi{\Km} \left[ \left(\PD_z\DM{u}\right)^2 + \left(\PD_z\DM{v}\right)^2 + \left(\PD_z\DM{w}\right)^2 \right] \\
\SDi{S_{\text{turb-transp}}}^{\phi\psi} & = - \PD_z (\rhoRef{} a_i \overline{w'_i\phi'_i\psi'_i}) \\
& = \PD_z (\rhoRef{} a_i \SDi{\Kh} \PD_z \IntraCVSDi{\phi}{\psi}) \\
\SDi{S_{\text{turb-transp}}}^{TKE} & = \PD_z (\rhoRef{} a_i \SDi{\Km} \PD_z \SDi{TKE}) \\
\SDi{S_{\text{dissip}}}
& = - \rhoRef{} a_i c_e \IntraCVSDi{\phi}{\psi} \frac{\SDi{TKE}^{1/2}}{\SDio{{l_{mix}}}}, \quad \text{Equation 38 in Tan et al.} \\
c_e & = 2 \\
\SDi{S_{\text{press}}}
& = - \aSDi{a} \left[ \IntraCVSDi{u}{(\partial_x p^{\dagger})} +
\IntraCVSDi{v}{(\partial_y p^{\dagger})} +
\IntraCVSDi{w}{(\partial_z p^{\dagger})}\right] \\
& = 0, \qquad \text{for now, need to derive correct formulation} \\
\SDi{S_{\text{buoyancy}}}^{TKE} & = \rhoRef{} \aSDi{a} \BuoyancyFlux \\
\SDi{S_{\text{MP-MSSP}}}^{\qt\qt}
& = \\
\SDi{S_{\text{MP-MSSP}}}^{\eint\eint}
& = \\
\end{align}\]
where additional variable definitions are in:
The following definitions are ordered in a dependency fashion; all variables are defined from variables already defined in previous subsections.
\[\begin{align}
c_K & = 0.1 \\
\text{tol}_{\InversionHeight\mathrm{-stable}} & = 0.01 \\
\end{align}\]
\[\begin{align}
\PhasePartition &= \{\qv, \ql, \qi\} \\
\qv &= \qt - \ql - \qi \\
\pvsat(T) &= \PTriple \left( \frac{T}{\TTriple} \right)^{\frac{\DeltaCp}{\Rv}} \exp{\frac{\RefLHv - \DeltaCp \TZero}{\Rv} \left( \frac{1}{\TTriple} - \frac{1}{T} \right)} \label{eq:pvsat} \\
\qvsat(T, \rho) &= \frac{\pvsat(T)}{\rho \Rv T} \label{eq:qvsat} \\
\qc &= \max(\qt - \qvsat, 0) \label{eq:qc} \\
\ql &= \lambda \qc \label{eq:ql} \\
\qi &= (1-\lambda) \qc \label{eq:qi} \\
\lambda(T) &= \Heaviside(T-\TFreeze) \label{eq:lambda} \\
\Heaviside &= \text{Heaviside function} \\
\end{align}\]
Functionally,
\[\begin{align}
\PhasePartition & = \PhasePartition(\qt, T, \rho) \\
\qvsat & = \qvsat(T, \rho) \\
\end{align}\]
\[\begin{align}
\EpsDV & = \frac{\Rv}{\Rd} \approx 1.61 \\
\EpsVD & = \frac{\Rd}{\Rv} \approx 0.62 \\
\Rm & = \Rd \left[1 + (\EpsDV-1) \qt - \EpsDV (\ql+\qi) \right] \\
\end{align}\]
\[\begin{align}
\Cvm &= (1 - \SDi{\qt}) \Cv{d} + \SDi{\qv} \Cv{v} + \SDi{\ql} \Cv{l} + \SDi{\qi} \Cv{i} \\
\Cpm &= (1 - \SDi{\qt}) \Cp{d} + \SDi{\qt} \Cp{v} \\
\end{align}\]
\[\begin{align}
\LatentHeatV{T} &= \RefLHv + (\Cp{v} - \Cp{l}) (T - \TTriple) \\
\LatentHeatS{T} &= \RefLHs + (\Cp{v} - \Cp{i}) (T - \TTriple) \\
\LatentHeatF{T} &= \RefLHf + (\Cp{l} - \Cp{i}) (T - \TTriple) \\
\end{align}\]
\[\begin{align}
\ExnerD(\pRef{}) = \left(\frac{\pRef{}}{\PTilde{}} \right)^{\Rd/\Cp{d}} \\
\ExnerM(\pRef{}, \PhasePartition) = \left(\frac{\pRef{}}{\PTilde{}} \right)^{\Rm/\Cpm} \\
\end{align}\]
where additional variable definitions are in:
Note that, while temperature may be computed using different thermodynamic formulations, ThermodynamicState's are immediately converted to the ($\qt, \eint, \rhoRef{}$)-formulation.
\[\begin{align}
\TDry & = \ThetaLiqIce \ExnerD \\
\end{align}\]
where additional variable definitions are in:
Functionally,
\[\begin{align}
\TDry & = \TDry(\ThetaLiqIce, \pRef) \\
\end{align}\]
Here, $T$ conditionally satisfies the non-linear set of equations, which can be solved using a standard root solver (e.g., Secant method):
\[\begin{align}
T =
\begin{cases}
\mathrm{satisfies} & \eint(T) = \Cvm (T - \TZero) + \qv \RefHintV - \qi \RefHintI & \qt > \qvsat(T, \rhoRef{}) \\
& \TZero + \frac{\eint(T)}{(1-\qt)\Cv{d} + \qt \Cv{v}} + \qt \RefHintV & \text{otherwise} \\
\end{cases}
\end{align}\]
where additional variable definitions are in:
Functionally,
\[\begin{align}
T & = T(\qt, \eint, \rhoRef) \\
\end{align}\]
Here, $T$ conditionally satisfies the non-linear set of equations, which can be solved using a standard root solver (e.g., Secant method):
\[\begin{align}
T =
\begin{cases}
\mathrm{satisfies} & \ThetaLiqIce \ExnerM = T \left(1 - \frac{ \RefLHv \ql + \RefLHs \qi}{\Cpm T} \right) & \qt > \qvsat(T, \rhoRef{}) \\
& \TDry & \text{otherwise} \\
\end{cases}
\end{align}\]
where additional variable definitions are in:
Functionally,
\[\begin{align}
T & = T(\qt, \ThetaLiqIce, \rhoRef, \pRef) \\
\end{align}\]
Using the hydrostatic balance, $\PD_z \pRef = - \rhoRef{} \grav$, and the ideal gas law, $\pRef = \rhoRef{} \Rm \TRef$, the reference state profiles are computed as:
\[\begin{align}
\PD_z \pRef & = - \grav \frac{\pRef}{\TRef \Rm} \\
\int_{\BC{\pRef}}^{\pRef} \frac{\TDry(\BC{\DM{\ThetaLiqIce}}, \pRef)}{\pRef} \PD \pRef & = - \frac{\grav}{\BC{\DM{\Rm}}} \int_{z_{min}}^{z} \PD z \\
\rhoRef{}(\pRef) & = \frac{\pRef{}}{\TDry(\BC{\DM{\ThetaLiqIce}}, \pRef) \BC{\DM{\Rm}}} \\
\alphaRef & = \frac{1}{\rhoRef{}(\pRef)} \\
\end{align}\]
where additional variable definitions are in:
\[\begin{align}\label{eq:MixingRatios}
r_{con} & = \frac{\qt+\ql}{1 - \qt} \\
r_{vap} & = \frac{\qt-\ql - \qi}{1 - \qt} \\
\end{align}\]
Fix: which virtual potential temperature is used
\[\begin{align}\label{eq:Theta}
\ThetaDry & = T/\ExnerD \\
\ThetaLiqIce & = \ThetaDry (1 - (\RefLHv \ql + \RefLHs \qi)/(\Cpm T)) \\
\ThetaVirt & = \ThetaDry (1 - r_{con} + 0.61 r_{vap}) \\
\ThetaVirt & = \theta \left(1 + 0.61 \qr - \ql \right) \\
\ThetaRho & = T \Rm/\ExnerD \\
\end{align}\]
where additional variable definitions are in:
\[\begin{align}\label{eq:ShearProduction}
|S|^2 &= (\PD_z \DM{u})^2 + (\PD_z \DM{v})^2 + (\PD_z \SDe{w})^2 \\
\end{align}\]
\[\begin{align}\label{eq:Buoyancy}
\SDi{b}^{\dagger} & = \grav (\SDi{\alpha} - \alphaRef)/\alphaRef \\
\SDi{b} & = \SDi{b}^{\dagger} - \sum_j a_j \SDj{b}^{\dagger} \\
\alpha_i & = \frac{\SDi{\Rm} \SDi{T}}{\pRef} \\
\end{align}\]
where additional variable definitions are in:
\[\begin{align}\label{eq:BuoyancyGradLong}
\SDi{\BuoyancyGrad} & = \PD_z \SDi{\ThetaLiqIce}
\left[ (1-f_c) \PD_{\ThetaLiqIce} b |_d + f_c \PD_{\ThetaLiqIce} b |_s \right] +
\PD_z \SDi{\qt} \left[ (1-f_c) \PD_{\qt} b |_d + f_c \PD_{\qt} b |_s \right] \\
f_c &= 0 \qquad \text{good for simple cases, need to confirm for more complex cases} \\
\PD_{\ThetaLiqIce} b |_d & = \frac{\grav}{\DM{\ThetaVirt}} \left[ 1 + \left( \frac{\Rv}{\Rd} - 1 \right) \SDi{\qt} \right] \\
\PD_{\ThetaLiqIce} b |_s &= \frac{\grav}{\DM{\ThetaVirt}} \left[ 1 + \frac{\Rv}{\Rd} \left(1 + \frac{\LatentHeatV{\SDi{T}}}{\Rv \SDi{T}} \right) \SDi{\qvsat} - \SDi{\qt} \right] \left( 1 + \frac{{\LatentHeatV{\SDi{T}}}^2}{\Cpm \Rv \SDi{T}^2} \SDi{\qvsat} \right)^{-1} \\
\PD_{\qt} b |_d &= \frac{\grav}{\DM{\ThetaVirt}} \left( \frac{\Rv}{\Rd} - 1 \right) \SDi{\ThetaDry} \\
\PD_{\qt} b |_s &= \left( \frac{{\LatentHeatV{\SDi{T}}}}{\Cpm \SDi{T}} \PD_{\ThetaLiqIce} b |_s - \frac{\grav}{\DM{\ThetaVirt}} \right) \SDi{\ThetaDry} \\
\end{align}\]
where additional variable definitions are in:
Pending.
Add definitions for universal functions (e.g., $\Psi_m$).
Variables in this section must be computed simultaneously because it requires the solution of a non-linear equation.
NOTE: All variables (Monin-Obhukov length, friction velocity, temperature scale) in Surface fluxes must be solved simultaneously
\[\begin{align}\label{eq:MOLen}
\MOLen = \begin{cases}
- \frac{\FrictionVelocity^3 \theta}{\VKConst \grav \SurfaceHeatFlux} & \SurfaceHeatFlux > 0 \\
0 & \text{otherwise} \\
\end{cases} \\
\end{align}\]
NOTE: All variables (Monin-Obhukov length, friction velocity, temperature scale) in Surface fluxes must be solved simultaneously
Knowns: $u_{\mathrm{ave}} = \sqrt{\DM{u}^2+\DM{v}^2}, \LayerThickness, \SurfaceRoughness{m}$
Unknowns: $\FrictionVelocity, \MOLen$, and $\SurfaceMomentumFlux$
\[\begin{align}\label{eq:FrictionVelocity}
u_{\mathrm{ave}} & = \frac{\FrictionVelocity}{\VKConst} \left[ \log\left(\frac{\LayerThickness}{\SurfaceRoughness{m}}\right) - \Psi_m\left(\frac{\LayerThickness}{\MOLen}\right) + \frac{\SurfaceRoughness{m}}{\LayerThickness} \Psi_m\left(\frac{\SurfaceRoughness{m}}{\MOLen}\right) + R_{z0m} \left\{ \psi_m\left(\frac{\SurfaceRoughness{m}}{\MOLen}\right) - 1 \right\} \right] \\
R_{z0m} & = 1 - \SurfaceRoughness{h}/\LayerThickness \\
\SurfaceMomentumFlux & = -\FrictionVelocity^2 , \label{eq:SurfaceMomentumFlux} \\
\end{align}\]
where $\Psi_m$ is defined in Appendix A, equations A6 in Nishizawa, S., and Y. Kitamura. "A Surface Flux Scheme Based on the Monin‐Obukhov Similarity for Finite Volume Models." Journal of Advances in Modeling Earth Systems 10.12 (2018): 3159-3175.
NOTE: All variables (Monin-Obhukov length, friction velocity, temperature scale) in Surface fluxes must be solved simultaneously
Knowns: $\theta_{\mathrm{ave}}, \theta_s, \LayerThickness, \SurfaceRoughness{h}$
Unknowns: $\FrictionVelocity, \MOLen$, and $\SurfaceHeatFlux$
\[\begin{align}\label{eq:TemperatureScale}
\theta_{\mathrm{ave}} - \theta_s & = \frac{Pr \TemperatureScale}{\VKConst} \left[ \log\left(\frac{\LayerThickness}{\SurfaceRoughness{h}}\right) - \Psi_h\left(\frac{\LayerThickness}{\MOLen}\right) + \frac{\SurfaceRoughness{h}}{\LayerThickness} \Psi_m\left(\frac{\SurfaceRoughness{h}}{\MOLen}\right) + R_{z0h} \left\{ \psi_h\left(\frac{\SurfaceRoughness{h}}{\MOLen}\right) - 1 \right\} \right] \\
R_{z0h} & = 1 - \SurfaceRoughness{h}/\LayerThickness \\
\SurfaceHeatFlux & = -\FrictionVelocity\TemperatureScale , \label{eq:SurfaceHeatFlux} \\
\end{align}\]
where $\Psi_h$ is defined in Appendix A, equation A6 in Nishizawa, S., and Y. Kitamura. "A Surface Flux Scheme Based on the Monin‐Obukhov Similarity for Finite Volume Models." Journal of Advances in Modeling Earth Systems 10.12 (2018): 3159-3175.
\[\begin{align}\label{eq:PrandtlNumber}
Pr_{neut} &= 0.74 \\
Pr(z) &= \begin{cases}
Pr_{neut} & \MOLen < 0 \\
Pr_{neut} \left[ \frac{1 + \omega_2 R_g - \sqrt{-4 R_g + (1+\omega_2 R_g)^2}}{2 R_g} \right] & \text{otherwise} \\
\end{cases} \\
\omega_2 &= \omega_1+1 \\
\omega_1 &= \frac{40}{13} \\
R_g &= \frac{\BuoyancyGrad}{|S|^2} \\
\end{align}\]
where additional variable definitions are in:
These mixing length have been tested for the environment, not the updrafts
\[\begin{align}\label{eq:MixingLength}
\SDio{{l_{mix}^m,}} &= \frac{\sum_j l_j e^{-l_j}}{\sum_j e^{-l_j}}, \qquad j = 1,2,3 \\
l_1 &= \frac{\sqrt{c_w\SDe{TKE}}}{\SDe{N}} \\
c_w &= 0.4 \\
\SDe{N} &= \frac{\grav \PD_z \SDe{\ThetaVirt}}{\SDe{\ThetaVirt}} , \qquad \text{(buoyancy frequency of environment)} \\
l_2 &= \frac{\VKConst z}{c_K \kappa^* \phi_m(z/\MOLen)} \\
\phi_m(\xi) &= \left( 1 + a_l \xi \right)^{-b_l} \\
(a_l, b_l) &=
\begin{cases}
(-100, 0.2) & \MOLen < 0 \\
(2.7, -1) & \text{otherwise} \\
\end{cases} \\
\kappa^* &= \frac{\FrictionVelocity}{\sqrt{\SDe{TKE}}} \\
l_3 &= \sqrt{\frac{c_{\varepsilon}}{c_K}} \sqrt{\SDe{TKE}}
\left[ \max(|S|^2 - \frac{1}{Pr(z)} \BuoyancyGrad, 0) \right]^{-1/2} \\
\end{align}\]
where additional variable definitions are in:
Smoothing function is provided in python file. The Prandtl number was used from Eq. 75 in Dan Li 2019 "Turbulent Prandtl number in the atmospheric BL - where are we now".
\[\begin{align}\label{eq:EddyDiffusivity}
\SDi{\Km} & = \begin{cases}
c_K \SDio{{l_{mix},}} \sqrt{\SDi{TKE}} & i = \iEnv \\
0 & \text{otherwise}
\end{cases} \\
\SDi{\Kh} & = \frac{\SDi{\Km}}{Pr} \\
\end{align}\]
where additional variable definitions are in:
Currently, $\BuoyancyFlux$ is hard-coded from the first expression (which was used in SCAMPy), however, this value should be computed from the SurfaceFluxes section.
\[\begin{align}\label{eq:BuoyancyFlux}
\SurfaceBuoyancyFlux & = \frac{\grav \BC{\alphaRef}}{\Cpm \BC{\SDi{T}}} (\SensibleSurfaceHeatFlux + (\EpsDV - 1) \Cpm \BC{\SDi{T}} \LatentSurfaceHeatFlux / \LatentHeatV{\BC{\SDi{T}}}) \\
\BuoyancyFlux & = - \SDi{\Kh} \SDi{\BuoyancyGrad} \\
\end{align}\]
Entrainment ($\epsilon_{i}$)
\[\begin{align}\label{eq:Entrainment}
\epsilon_{i} &= c_{\epsilon} \frac{\max(\SDi{b}, 0)}{\SDi{w}^2} \\
c_{\epsilon} &= 0.12 \\
\end{align}\]
Detrainment ($\delta_{j}$):
\[\begin{align}\label{eq:Detrainment}
\delta_{i} &= c_{\delta} \frac{|\min(\SDi{b}, 0)|}{\SDi{w}^2} + \delta_{B} \Heaviside(\SDi{\ql}) \\
c_{\delta} &= c_{\delta,0} + \Gamma(\aSDi{a}) \\
\Gamma(\aSDi{a}) &= 0 \\
c_{\delta,0} &= c_{\epsilon} = 0.12 \\
\delta_B &= 0.004 [m^{-1}] \\
\Heaviside &= \text{Heaviside function} \\
\end{align}\]
where additional variable definitions are in:
\[\begin{align}\label{eq:InversionHeight}
\SDio{\InversionHeight} &=
\begin{cases}
\left[ (\PD_z \ThetaRho)^{-1} (\BC{\ThetaRho} - \ThetaRho|_{z_1}) + z_1 \right] & \simparam{\BC{\DM{u}}}^2 + \simparam{\BC{\DM{v}}}^2 <= \text{tol}_{\InversionHeight\mathrm{-stable}} \\
\left[ (\PD_z Ri_{bulk})^{-1} (\hyperparam{Ri_{bulk, crit}} - Ri_{bulk}|_{z_2}) + z_2 \right] & \text{otherwise} \\
\end{cases} \\
z_1 &= \min_z (\ThetaRho(z) > \BC{\ThetaRho}) \\
z_2 &= \min_z (Ri_{bulk}(z) > \hyperparam{Ri_{bulk, crit}}) \\
Ri_{bulk} &= \grav z \frac{(\ThetaRho/\BC{\ThetaRho} - 1)}{\simparam{\DM{u}}^2 + \simparam{\DM{v}}^2} \\
\end{align}\]
where additional variable definitions are in:
\[\begin{align}\label{eq:ConvectiveVelocity}
\SDio{\ConvectiveVelocity} &= (\max(\BuoyancyFlux \SDio{\InversionHeight}, 0))^{1/3} \\
\end{align}\]
where additional variable definitions are in:
\[\begin{align}\label{eq:MixingLengthOld}
\SDio{{l_{mix},}} &= (l_A^{-1} + l_B^{-1})^{-1} \\
l_A &= \VKConst z \left( 1 + a_l \frac{z}{\MOLen} \right)^{b_l} \\
\SDio{{l_B}} &= \SDio{\tau} \SDi{TKE} \\
(a_l, b_l) &=
\begin{cases}
(-100, 0.2) & \MOLen < 0 \\
(2.7, -1) & \text{otherwise} \\
\end{cases} \\
\SDio{\tau} &= \SDio{\InversionHeight}/\SDio{\ConvectiveVelocity} \\
\end{align}\]
where additional variable definitions are in:
Here, we specify boundary conditions (BCs) by their type, Dirichlet (D) or Neumann (N), and their value.
\[\begin{align}
\Gamma_{\phi}(F_1, F_2)
& = \begin{cases}
4 \frac{F_1 F_2}{\FrictionVelocity^2} (1 - 8.3\zLL/\MOLen)^{-2/3} & \MOLen < 0 \\
4 \frac{F_1 F_2}{\FrictionVelocity^2} & \text{otherwise}
\end{cases} \\
\Gamma_{TKE}
& = \begin{cases}
3.75 {\FrictionVelocity}^2 + 0.2 {\ConvectiveVelocity}^2 + {\FrictionVelocity}^2 (-\zLL/\MOLen)^{2/3} & \MOLen < 0 \\
3.75 {\FrictionVelocity}^2 & \text{otherwise}
\end{cases} \\
\SensibleSurfaceHeatFlux & = \BC{\TCV{w}{\eint}} \Cpm \rhoRef \\
\LatentSurfaceHeatFlux & = \BC{\TCV{w}{\qt}} \LatentHeatV{T} \rhoRef \\
F_{\eint}(\SensibleSurfaceHeatFlux) & = \frac{\SensibleSurfaceHeatFlux}{\Cpm} = \BC{\TCV{w}{\eint}} \rhoRef \\
F_{\qt}(\LatentSurfaceHeatFlux) & = \frac{\LatentSurfaceHeatFlux}{\LatentHeatV{T}} = \BC{\TCV{w}{\qt}} \rhoRef \\
\end{align}\]
where additional variable definitions are in:
and equation \eqref{eq:TopPercentile} represents the mean of the top $x$-fraction of a standard normal distribution (Neggers et al., 2009).
\[\begin{align}
\Phi^{-1}(x) &= \text{inverse cumulative distribution function}, \label{eq:InverseCDF} \\
\mathcal D(x) &= \frac{1}{\sqrt{2\pi x}} \exp{- \frac{1}{2} (\Phi^{-1}(1-x))^2 } , \label{eq:TopPercentile} \\
\end{align}\]
\[\begin{align}
c_{frac} = 0.1, \quad
\BCB{\aSDi{a}} =
\begin{cases}
1-c_{frac}, & i = \iEnv{} \\
\frac{c_{frac}}{\Nsd}, & i \ne \iEnv{}
\end{cases}, \quad
\BCT{\aSDi{a}} =
\begin{cases}
1-c_{frac}, & i = \iEnv{} \\
\frac{c_{frac}}{\Nsd}, & i \ne \iEnv{}
\end{cases}
\end{align}\]
Top boundary
\[\begin{align}
\BCT{\SDi{w}} &= 0 \\
\PD_z \BCT{\SDi{\qt}} &= 0 \\
\PD_z \BCT{\SDi{\eint}} &= 0 \\
\end{align}\]
Bottom boundary
Need value for $C_{\eint}$.
\[\begin{align}
\BCB{\SDi{w}} &= 0 \\
- \SDi{\Kh} \PD_z \BCB{\SDi{\qt}} &= \TCV{w}{\qt} + \mathcal D(\aSDi{a}) \sqrt{C_{\qt}^2 \WindSpeed^2\Gamma_{\phi}(\TCV{w}{\qt} , \TCV{w}{\qt} )} \\
- \SDi{\Kh} \PD_z \BCB{\SDi{\eint}} &= \TCV{w}{\eint} + \mathcal D(\aSDi{a}) \sqrt{C_{\eint}^2 \WindSpeed^2\Gamma_{\phi}(\TCV{w}{\eint}, \TCV{w}{\eint} )} \\
C_{\qt} &= 0.001133 \\
C_{\ThetaLiqIce} &= 0.001094 \\
C_{\eint} &= \\
\end{align}\]
where additional variable/function definitions are in:
Top boundary
\[\begin{align}
\BCT{\SDi{TKE}} & = 0 \\
\PD_z \BCT{\IntraCVSDi{\qt}{\qt}} & = 0 \\
\PD_z \BCT{\IntraCVSDi{\eint}{\eint}} & = 0 \\
\PD_z \BCT{\IntraCVSDi{\eint}{\qt}} & = 0 \\
\end{align}\]
Currently, we only account for the intra sub-domain covariance, but we would like to also account for the inter sub-domain covariance for all but the $TKE$.
Bottom boundary
\[\begin{align}
\BCB{\SDi{TKE}} & = \Gamma_{TKE} \\
\BCB{\IntraCVSDi{\qt}{\qt}} & = \Gamma_{\phi}(\TCV{w}{\qt} , \TCV{w}{\qt} ) \\
\BCB{\IntraCVSDi{\eint}{\eint}} & = \Gamma_{\phi}(\TCV{w}{\qt} , \TCV{w}{\eint} ) \\
\BCB{\IntraCVSDi{\eint}{\qt}} & = \Gamma_{\phi}(\TCV{w}{\eint}, \TCV{w}{\eint} ) \\
\end{align}\]
where additional variable/function definitions are in:
- BC functions $\Gamma_{TKE}$, $\Gamma_{\phi}$, $F_{\eint}$, $\SensibleSurfaceHeatFlux$, $F_{\qt}$, $\LatentSurfaceHeatFlux$.