brodiepearson/mpas-model

Repository for MPAS models and shared framework releases.

mpas-model's Issues

Comment & rearrange tendency budgets in ADC closure

The following ADC tendency budgets require:

• Commenting to make the meaning of each term clearer
• Reorganization so that the individual tendency terms diagnosed by the ADC scheme are identical to the tendency terms that we get calculate in LES or the standard RANS budgets.

For example with reogranization, sometimes the parameterized pressure terms are spread across multiple ADC tendency terms where they are mixed with other processes, or they appear implicitly (in the case of the slow Rotta term). In contrast, LES and RANS they are associated with a single pressure-strain term in each budget.

The budgets that need to be commented/reorganized are

• w2
• wt
• ws
• w3
• u2
• v2
• uw
• vw
• uv
• ut
• vt
• us
• vs

After this, perhaps we should also comment the sub-plume scale budgets:

• KspsUtend
• KspsDtend
• epstend

Test issue

I think there are a couple of issues in the sub-plume terms in the $\overline{w'w'}$, $\overline{w'T'}$ and $\overline{w'S'}$ budgets, but I first wanted to double check that these are typos rather than me simply missing some other element of the MPAS scheme that makes these formulations correct. Changing these terms doesn't have a significant impact in the convective cases, but it might be important when the sub-plume terms play a more active role (wind-driven mixing?).

The Lappen and Randall papers present the sub-plume terms for the $\overline{w'w'}$ budget (Eqs. 40 and 48 of their first paper) and the $\overline{w'T'}$ budget (Eq. 49 of their first paper and Eq. 27 of their second paper). Noting that our ocean ADC scheme uses $\sigma$ to represent downwelling plume area rather than the upwelling plume area used in LR papers, these equations in the ocean ADC scheme are:

$$\left(\frac{\partial \overline{w'w'}}{\partial t}\right)^{sps} = -2M_c\left[ \frac{1}{1-\sigma}\frac{\partial}{\partial z} \left([1-\sigma] \overline{w'w'}^{sps,up}\right) - \frac{1}{\sigma}\frac{\partial}{\partial z} \left(\sigma \overline{w'w'}^{sps,down}\right) \right] $$

$$\left(\frac{\partial \overline{w'T'}}{\partial t}\right)^{sps} = \left(\frac{\partial \overline{w'T'}}{\partial t}\right)^{sps, w} + \left(\frac{\partial \overline{w'T'}}{\partial t}\right)^{sps, T}$$

where

$$\left(\frac{\partial \overline{w'T'}}{\partial t}\right)^{sps, w} = -\sigma(1-\sigma)(T_u-T_d)\left[ \frac{1}{1-\sigma}\frac{\partial}{\partial z} \left([1-\sigma] \overline{w'w'}^{sps,up}\right) - \frac{1}{\sigma}\frac{\partial}{\partial z} \left(\sigma \overline{w'w'}^{sps,down}\right) \right] $$

$$\left(\frac{\partial \overline{w'T'}}{\partial t}\right)^{sps, T} = \frac{M_c}{1-\sigma}\frac{\partial}{\partial z} \left(K_h[1-\sigma] \frac{\partial \overline{T}}{\partial z}\right)^{up} - \frac{M_c}{\sigma}\frac{\partial}{\partial z} \left(K_h\sigma \overline{T}{\partial z}\right)^{down} $$

Typos in the ADC scheme?

The ADC scheme looks like it is applying sub-plume terms that are different from the above formulations. The $\overline{w'w'}$ term seems to have the wrong sign, while the $\overline{w'T'}$ and $\overline{w's'}$ terms look like the LR formulation where $\sigma$ is the upwelling rather than downwelling plume area. The ocean ADC formulations are (with differences highlighted in red and the associated ADC code):

$$\left(\frac{\partial \overline{w'w'}}{\partial t}\right)^{sps} = \color{red}+2M_c\left[ \frac{1}{1-\sigma}\frac{\partial}{\partial z} \left([1-\sigma] \overline{w'w'}^{sps,up}\right) - \frac{1}{\sigma}\frac{\partial}{\partial z} \left(\sigma \overline{w'w'}^{sps,down}\right) \right] $$

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Line 603 in aff0c19

w2tend5(k,iCell) = Mc(k,iCell)*(Swumd(k-1,iCell) + Swumd(k,iCell))

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Lines 459 to 462 in aff0c19

	Swumd(k,iCell) = - 2.0_RKIND/3.0_RKIND*(1.0_RKIND/sigav*(areaFraction(k,iCell)*KspsD(i1,k,iCell) - &
	areaFraction(k+1,iCell)*KspsD(i1,k+1,iCell)) / dz - 1.0_RKIND / (1.0_RKIND - &
	sigav)*((1.0_RKIND - areaFraction(k,iCell))*KspsU(i1,k,iCell) - (1.0_RKIND - &
	areaFraction(k+1,iCell))*KspsU(i1,k+1,iCell)) / dz)

$$\left(\frac{\partial \overline{w'T'}}{\partial t}\right)^{sps, w} = \color{red}+\sigma(1-\sigma)(T_u-T_d)\left[ \frac{1}{1-\sigma}\frac{\partial}{\partial z} \left([1-\sigma] \overline{w'w'}^{sps,up}\right) - \frac{1}{\sigma}\frac{\partial}{\partial z} \left(\sigma \overline{w'w'}^{sps,down}\right) \right] $$

$$\left(\frac{\partial \overline{w'T'}}{\partial t}\right)^{sps,T} = \frac{M_c}{\color{red}{\sigma}}\frac{\partial}{\partial z} \left(K_h\color{red}{\sigma} \frac{\partial \overline{T}}{\partial z}\right)^{up} - \frac{M_c}{\color{red}{1-\sigma}}\frac{\partial}{\partial z} \left(K_h\color{red}{(1-\sigma)} \frac{\partial \overline{T}}{\partial z}\right)^{down} $$

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Lines 652 to 658 in aff0c19

	wttend6(k,iCell) = areaFraction(k,iCell)*(1.0_RKIND - areaFraction(k,iCell))* &
	tumd(k,iCell)*0.5_RKIND*(Swumd(k-1,iCell) + Swumd(k,iCell)) - Mc(k,iCell)* &
	(1.0_RKIND / areaFraction(k,iCell) * (areaFraction(k-1,iCell)*wt_spsU(k-1,iCell) - &
	areaFraction(k+1,iCell)*wt_spsU(k+1,iCell)) / (ze(k-1,iCell) - ze(k+1,iCell)) &
	-1.0_RKIND/(1.0_RKIND - areaFraction(k,iCell))*((1.0_RKIND - areaFraction(k-1,iCell))* &
	wt_spsD(k-1,iCell) - (1.0_RKIND - areaFraction(k+1,iCell))*wt_spsD(k+1,iCell)) / &
	(ze(k-1,iCell) - ze(k+1,iCell)))

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Line 1038 in aff0c19

wt_spsU(k,iCell) = -KhU(k,iCell)*Tz

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Line 1043 in aff0c19

wt_spsD(k,iCell) = -KhD(k,iCell)*Tz

Fixes

If the above is correct, the easiest fix would be to

1. Change the sign of Swumd and
2. Flip $\sigma$ and $(1-\sigma)$ usage in the $\left(\frac{\partial \overline{w'T'}}{\partial t}\right)^{sps,T}$ code of wttend6 and wstend6

Some typos in sub-plume and entrainment terms

I think there are a couple of typos in several terms that need to be fixed. These typos appeared because of the switch from atmosphere to ocean orientation, which has led to some terms being flipped correctly, and others not.

Typos I potentially identified:

• The $\overline{w'\theta'}$ sub-plume scale term incorrectly associates areaFraction (downwelling area) with wt_spsU (upwelling flux), and vice versa for downwelling flux. In the terms containing wt_spsU and wt_spsD you should swap around $\sigma$ and $(1-\sigma)$ so that $\sigma$ is associated with wt_spsD and vice versa.
  

  MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

  Lines 652 to 658 in 1a3f760

  	wttend6(k,iCell) = areaFraction(k,iCell)*(1.0_RKIND - areaFraction(k,iCell))* &
  	tumd(k,iCell)*0.5_RKIND*(Swumd(k-1,iCell) + Swumd(k,iCell)) - Mc(k,iCell)* &
  	(1.0_RKIND / areaFraction(k,iCell) * (areaFraction(k-1,iCell)*wt_spsU(k-1,iCell) - &
  	areaFraction(k+1,iCell)*wt_spsU(k+1,iCell)) / (ze(k-1,iCell) - ze(k+1,iCell)) &
  	-1.0_RKIND/(1.0_RKIND - areaFraction(k,iCell))*((1.0_RKIND - areaFraction(k-1,iCell))* &
  	wt_spsD(k-1,iCell) - (1.0_RKIND - areaFraction(k+1,iCell))*wt_spsD(k+1,iCell)) / &
  	(ze(k-1,iCell) - ze(k+1,iCell)))

• Same as above, but with $\overline{w's'}$.
  

  MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

  Lines 708 to 714 in 1a3f760

  	wstend6(k,iCell) = areaFraction(k,iCell)*(1.0_RKIND - areaFraction(k,iCell))* &
  	sumd(k,iCell)*0.5_RKIND*(Swumd(k-1,iCell) + Swumd(k,iCell)) - Mc(k,iCell)* &
  	(1.0_RKIND / areaFraction(k,iCell) * (areaFraction(k-1,iCell)*ws_spsU(k-1,iCell) - &
  	areaFraction(k+1,iCell)*ws_spsU(k+1,iCell)) / (ze(k-1,iCell) - ze(k+1,iCell)) &
  	-1.0_RKIND/(1.0_RKIND - areaFraction(k,iCell))*((1.0_RKIND - areaFraction(k-1,iCell))* &
  	ws_spsD(k-1,iCell) - (1.0_RKIND - areaFraction(k+1,iCell))*ws_spsD(k+1,iCell)) / &
  	(ze(k-1,iCell) - ze(k+1,iCell)))

• I think there is an issue with the sign of Swumd, where it should be upwelling minus downwelling, but it has been defined the other way around. Remove the minus sign at the begining of this line:
  

  MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

  Lines 459 to 462 in 1a3f760

  	Swumd(k,iCell) = - 2.0_RKIND/3.0_RKIND*(1.0_RKIND/sigav*(areaFraction(k,iCell)*KspsD(i1,k,iCell) - &
  	areaFraction(k+1,iCell)*KspsD(i1,k+1,iCell)) / dz - 1.0_RKIND / (1.0_RKIND - &
  	sigav)*((1.0_RKIND - areaFraction(k,iCell))*KspsU(i1,k,iCell) - (1.0_RKIND - &
  	areaFraction(k+1,iCell))*KspsU(i1,k+1,iCell)) / dz)

• There may be a sign error in the $\overline{w'^3}$ entrainment/detrainment terms (I don't think there is, but my calculations got confusing so it would be nice to check). You should check whether w3tend1 has the same sign as $\overline{w'^3}$ or the opposite sign (it should be opposite!). If it w3tend1 and $\overline{w'^3}$ have the same sign, add a minus sign at the beginning of the following line.
  

  MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

  Lines 490 to 491 in 1a3f760

  	w3tend1(k,iCell) = -wumdav**3.0_RKIND*(Eav*(3.0_RKIND*sigav - 2.0_RKIND) &
  	+ Dav*(3.0_RKIND*sigav - 1.0_RKIND))

Physical Reasoning

I think these issues stem from conversion from atmosphere to ocean systems (intense upwelling to intense downwelling). Particularly, areaFraction in the new scheme represents downwelling area as it is <0.5 when skewness is negative due to ocean convection, versus upwelling area in the Lappen and Randall papers.

MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

Line 364 in 1a3f760

lareaFraction = 0.5_RKIND + 0.5_RKIND*Sw / sqrt(4.0_RKIND + Sw**2)

The one place where the choice of $\sigma$ for downwelling does affect results is:

1. in a few of the sub-plume terms where the area must be coupled with the correct plume variable in the derivatives. This is done correctly for the TKE based sub-plume terms, but is not correct for the tracer flux subplume terms.
   

   MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

   Lines 459 to 462 in 1a3f760

   	Swumd(k,iCell) = - 2.0_RKIND/3.0_RKIND*(1.0_RKIND/sigav*(areaFraction(k,iCell)*KspsD(i1,k,iCell) - &
   	areaFraction(k+1,iCell)*KspsD(i1,k+1,iCell)) / dz - 1.0_RKIND / (1.0_RKIND - &
   	sigav)*((1.0_RKIND - areaFraction(k,iCell))*KspsU(i1,k,iCell) - (1.0_RKIND - &
   	areaFraction(k+1,iCell))*KspsU(i1,k+1,iCell)) / dz)

   
   vs.
   

   MPAS-Model/src/core_ocean/shared/mpas_ocn_adcReconstruct.F

   Lines 653 to 658 in 1a3f760

   	tumd(k,iCell)*0.5_RKIND*(Swumd(k-1,iCell) + Swumd(k,iCell)) - Mc(k,iCell)* &
   	(1.0_RKIND / areaFraction(k,iCell) * (areaFraction(k-1,iCell)*wt_spsU(k-1,iCell) - &
   	areaFraction(k+1,iCell)*wt_spsU(k+1,iCell)) / (ze(k-1,iCell) - ze(k+1,iCell)) &
   	-1.0_RKIND/(1.0_RKIND - areaFraction(k,iCell))*((1.0_RKIND - areaFraction(k-1,iCell))* &
   	wt_spsD(k-1,iCell) - (1.0_RKIND - areaFraction(k+1,iCell))*wt_spsD(k+1,iCell)) / &
   	(ze(k-1,iCell) - ze(k+1,iCell)))

Add Stokes terms into budgets

The effects of Stokes drift on upper ocean mixing need to be incorporated in the ADC scheme. These contributions can be separated into

1. Explicit forces based on the Craik-Leibovich (CL/wave-averaged) equations
2. Closure terms modified by the Stokes drift

Explicit Forces

The wave-averaged equations introduce additional forces due to the Stokes drift ($\mathbf{u^{st}}$) the Stokes vortex force, the Stokes-Coriolis force, and a Bernouilli pressure modification due to Stokes drift. Together these terms modify the momentum equations:

$\frac{\partial u_i}{\partial t}+u_k\frac{\partial u_i}{\partial x_k}=\ldots - \frac{1}{2}\frac{\partial (u_k^{st})^2}{\partial x_i} - \frac{1}{2}\frac{\partial (u_k^{st}u_k)}{\partial x_i}+ \epsilon_{ijk}\epsilon_{klm}u_j^{st}\frac{\partial u_m}{\partial x_l} -\epsilon_{ijk}f_ju_k^{st}$.

We are interested in how these forces affect the turbulence statistics (fluxes, variances etc.). The quadratic Stokes term and the Coriolis-Stokes term disappear in the fluctuating momentum budgets (assuming the Stokes term is steady in time and horizontal space). The only effect that Stokes drift has in tracer ($\phi$) budgets is through the advection of that tracer.

Closure Terms

The pressure field within a turbulent flow is coupled to the momentum equations. By taking the divergence of the (vector) momentum budget and using incompressibility, an equation relating the turbulent pressure field to the other forces that drive the turbulent flow can be derived. This means that each force in the flow has an imprint on the turbulent pressure field. As a result the pressure-strain ($\Pi_{ij} = -\overline{u_i'[\partial p'/ \partial x_j]} - \overline{u_j'[\partial p'/ \partial x_i]}$) and pressure-scalar ($\Pi_{i\phi}=- \overline{\phi'[\partial p'/ \partial x_i]}$) terms that appear in Reynolds-averaged turbulence budgets will be affected by all the forces.

Pearson et al., 2019 parameterized the effect that Stokes drift has on the pressure-strain terms in Langmuir turbulence, using a suite of LES experiments. This parameterization was:
$\Pi_{ij}^{Stokes} = \alpha_0^{st}e S_{ij}^{st} + \alpha_1^{st}\left(S_{ik}^{st}a_{kj}+S_{jk}^{st}a_{ki}-(2/3)\delta_{ij}S_{kl}^{st}a_{kl}\right) + \alpha_2^{st}\left(W_{ik}^{st}a_{kj}+W_{jk}^{st}a_{ki}\right),$

where $S_{ij}^{st}=\frac{1}{2}\left[\frac{\partial u_{i}^{st}}{\partial x_j} +\frac{\partial u_{j}^{st}}{\partial x_i}\right]$, and $W_{ij}^{st}=\frac{1}{2}\left[\frac{\partial u_{i}^{st}}{\partial x_j} -\frac{\partial u_{j}^{st}}{\partial x_i}\right]$ are the strain rate and vorticity respectively. The $\alpha^{st}$ terms are constants, and values for these constants are suggested in the paper.

The pressure-scalar terms have not been investigated in wave-affected flows. The only attempt to parameterize them was in Harcourt, 2013 where they were assumed to be identical in form to the Eulerian shear-based pressure-scalar terms:

$$\Pi_{i \phi}^{rapid} = \begin{pmatrix} \alpha_{\phi} \overline{w' \phi'}\frac{\partial u}{\partial z} \\\ \alpha_{\phi} \overline{w' \phi'}\frac{\partial v}{\partial z} \\\ 0 \end{pmatrix} \quad \textrm{and} \quad \Pi_{i \phi}^{stokes} = \begin{pmatrix} \alpha_{\phi}^{st} \overline{w' \phi'}\frac{\partial u^{st}}{\partial z} \\\ \alpha_{\phi}^{st} \overline{w' \phi'}\frac{\partial v^{st}}{\partial z} \\\ 0 \end{pmatrix}$$

I suspect that the above pressure-scalar closure is not appropriate for Langmuir turbulence, as it simply replicates exactly the same functional form for both the shear-driven and Stokes-driven pressure effects, with Eulerian and Stokes velocities being interchanged between the two.

However a closer look at the the tracer flux budgets (below) shows that the rapid closure exactly cancels a fraction of the shear-production of tracer fluxes. There is no shear production of the vertical tracer flux ($\overline{w'\phi'}$) but there can be significant production of horizontal fluxes ($\overline{u'\phi'}$ and $\overline{v'\phi'}$). In contrast, Stokes drift produces vertical fluxes of tracers, but it does not produce horizontal tracer fluxes. A more appropriate pressure-scalar closure for Stokes drift effects, which cancels this Stokes production term, would be,

$$\Pi_{i \phi}^{rapid} = \begin{pmatrix} \alpha_{\phi} \overline{w' \phi'}\frac{\partial u}{\partial z} \\\ \alpha_{\phi} \overline{w' \phi'}\frac{\partial v}{\partial z} \\\ 0 \end{pmatrix} \quad \textrm{and} \quad \Pi_{i \phi}^{stokes} = \begin{pmatrix} 0 \\\ 0 \\\ \alpha_{\phi}^{st} \left(\overline{u' \phi'}\frac{\partial u^{st}}{\partial z} + \overline{v' \phi'}\frac{\partial v^{st}}{\partial z}\right) \end{pmatrix}$$

We will use the simplest form of the pressure-scalar closure that allows both of these closure options, which is based on Mironov, 2001 . This closure requires two closure constants for Stokes pressure-scalar closure (and analogous constants for the rapid closure, which is already coded into the ADC scheme):

$$ \Pi_{i \phi}^{stokes} = \begin{pmatrix} \left(\frac{\alpha_{\phi 1}^{st} + \alpha_{\phi 2}^{st}}{2}\right) \overline{w' \phi'}\frac{\partial \overline{u^{st}}}{\partial z} \\\ \left(\frac{\alpha_{\phi 1}^{st} + \alpha_{\phi 2}^{st}}{2}\right) \overline{w' \phi'}\frac{\partial \overline{v^{st}}}{\partial z} \\\ \left(\frac{\alpha_{\phi 1}^{st} - \alpha_{\phi 2}^{st}}{2}\right) \left(\overline{u' \phi'}\frac{\partial u^{st}}{\partial z} + \overline{v' \phi'}\frac{\partial v^{st}}{\partial z}\right) \end{pmatrix}$$

where the Harcourt 2013 closure would correspond to $\alpha_{\phi 1}=\alpha_{\phi 2}=\alpha/2$ and the production-cancellation option would correspond to $\alpha_{\phi 1}=-\alpha_{\phi 2}=\alpha/2$.

This form of closure for the Stokes pressure-scalar term also allows a range of other closures. For example, Canuto et al. (2001) suggested that flow constraints and idealized turbulence experiments suggest that the rapid pressure-scalar term should be parameterized as above but with $\frac{4}{3} \alpha_{\phi 1} = \frac{4}{5} \alpha_{\phi 2}$. An analogous form for the Stokes term is allowed by the above closure.

It is not clear how the Stokes drift could affect the pressure term in the skewness ($\overline{w'w'w'}$) budget, so we will neglect it for now.

Summary of Stokes drift effects in ADC equations

The equations governing the turbulence statistics of a flow can be found by Reynolds averaging to find equations for turbulent statistics. The inclusion of a Stokes drift, both explicitly and in pressure-related closure assumptions, leads to the following additional terms in the budgets for turbulent fluxes and variances (note that $D/Dt$ includes the advection by both Eulerian flow and Stokes drift and bold terms denote the explicit (non-closure) Stokes terms),

Variances & Skewness

$\frac{D_L \overline{w'w'}}{D t} =\ldots - \left(\mathbf{2} - \frac{\alpha^{st}_1}{3}+\alpha_2^{st} \right)\left[\overline{u'w'}\frac{\partial u^{st}}{\partial z} + \overline{v'w'}\frac{\partial v^{st}}{\partial z}\right],$

$\frac{D_L \overline{u'u'}}{D t} =\ldots + \left(\frac{\alpha_1^{st}}{3} + \alpha_2^{st}\right)\overline{u'w'}\frac{\partial u^{st}}{\partial z} - \frac{2\alpha_1^{st}}{3}\overline{v'w'}\frac{\partial v^{st}}{\partial z}$

$\frac{D_L \overline{v'v'}}{D t} =\ldots + \left(\frac{\alpha_1^{st}}{3} + \alpha_2^{st}\right)\overline{v'w'}\frac{\partial v^{st}}{\partial z} - \frac{2\alpha_1^{st}}{3}\overline{u'w'}\frac{\partial u^{st}}{\partial z}$

$\frac{D_L \overline{\phi'\phi'}}{D t} =\ldots \textrm{No wave or pressure effects}$

$\frac{D_L \overline{w'w'w'}}{D t} =\ldots - 3\overline{u'w'w'}\frac{\partial u^{st}}{\partial z} - 3\overline{v'w'w'}\frac{\partial v^{st}}{\partial z}, \qquad \textrm{ignoring Stokes effects in $- 3\overline{w'w'\frac{\partial p'}{\partial z}}$}$

Momentum Fluxes

$\frac{D_L \overline{u'w'}}{D t} =\ldots - \left(\mathbf{1} - \frac{\alpha_1^{st} - \alpha_2^{st}}{2}\right)\left[\overline{u'u'}\frac{\partial u^{st}}{\partial z} + \overline{u'v'}\frac{\partial v^{st}}{\partial z}\right] +\left(\frac{\alpha_0^{st}}{2} - \frac{2\alpha_1^{st}}{3}\right)e\frac{\partial u^{st}}{\partial z} + \left(\frac{\alpha_1^{st} + \alpha_2^{st}}{2}\right)\overline{w'w'}\frac{\partial u^{st}}{\partial z}$

$\frac{D_L \overline{v'w'}}{D t} =\ldots - \left(\mathbf{1} - \frac{\alpha_1^{st} - \alpha_2^{st}}{2}\right)\left[\overline{u'v'}\frac{\partial u^{st}}{\partial z} + \overline{v'v'}\frac{\partial v^{st}}{\partial z}\right] + \left(\frac{\alpha_0^{st}}{2} - \frac{2\alpha_1^{st}}{3}\right)e\frac{\partial v^{st}}{\partial z} + \left(\frac{\alpha_1^{st} + \alpha_2^{st}}{2}\right)\overline{w'w'}\frac{\partial v^{st}}{\partial z}$

$\frac{D_L \overline{u'v'}}{D t} =\ldots + \frac{\alpha_1^{st}+\alpha_2^{st}}{2} \left(\overline{v'w'}\frac{\partial u^{st}}{\partial z} + \overline{u'w'}\frac{\partial v^{st}}{\partial z} \right)$

Tracer fluxes

$ \frac{D_L \overline{w'\phi'}}{D t} =\ldots - \left[\mathbf{1}-\left(\frac{\alpha_{\phi 1}^{st} - \alpha_{\phi 2}^{st}}{2}\right)\right]\left[\overline{u'\phi'}\frac{\partial u^{st}}{\partial z} + \overline{v'\phi'}\frac{\partial v^{st}}{\partial z}\right]$

$\frac{D_L \overline{u'\phi'}}{D t} =\ldots + \left(\frac{\alpha_{\phi 1}^{st} + \alpha_{\phi 2}^{st}}{2}\right)\overline{w'\phi'}\frac{\partial\overline{u^{st}}}{\partial z}$

$\frac{D_L \overline{v'\phi'}}{D t} =\ldots + \left(\frac{\alpha_{\phi 1}^{st} + \alpha_{\phi 2}^{st}}{2}\right)\overline{w'\phi'}\frac{\partial\overline{v^{st}}}{\partial z}$