Typos in sub-plume scale terms

[BrodiePearson]

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

[vanroekel]

@BrodiePearson any chance you have a latex version of this? I'm struggling to get mathjax going.

[vanroekel]

my first gut feeling though is that I don't doubt this is wrong, but am pretty doubtful it will matter, even for wind driven cases, but could be wrong for sure

[vanroekel]

@BrodiePearson I got the mathjax to work. so disregard above. The one critical assumption we need to verify if the code currently does consider \sigma to be the downwelling plume. I'm not sure that is true. I don't know if we actually switched around. My recollection is we never switched. Is this something you all can check?

[BrodiePearson]

@vanroekel It looks like it was switched round to be downwelling plume area in this PR in 2021.

The Lappen and Randall formulation for upwelling plume area is:

The ocean ADC code is

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

Line 364 in a8d1ec4

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

The ocean ADC code flips the middle sign to a plus, which makes lareaFraction the downwelling area fraction. This makes conceptual sense because it makes lareaFraction small (<0.5) when skewness is negative (strong downwelling plumes) and large (>0.5) when skewness is positive (strong upwelling plumes with weak broad downwelling)

[vanroekel]

Thanks for checking into it @BrodiePearson, given this I fully agree with your initial comment that these are indeed bugs. Are you or Amrapalli able to test the fixes? I'm not sure when Filipe or I can get to them

[amrapallig]

@BrodiePearson and @vanroekel
I checked by fixing the suggested bugs for w2tend5 and wttend6 terms. Here are two plots for convective (upper panel) and wind case (lower panel) for 2days output. It seems the wind case is slightly affected and convection case shows no significant changes. Again wind case still needs tuning.

Few points to note here:

1. For our draft, I had used old commit (287e603) before adding Stokes terms and following ADC results are with latest commit (aff0c19). In addition to Stokes terms the latest commit also fixed two bugs in time stepping and uwtend, vwtend terms.

2. Blue and cyan are with the same namelist parameters as in the draft and those show over-mixing with updated commit. Hence I tuned CwwE and CwwD parameters to 1 and 1.5 which compare best to both NCAR and Oceananigan LES (red and magenta). all ADC cases need to use updated parameters.

3. Red dashed lines are with updated parameters and without above suggested sub-plume scale corrections. Magenta dotted lines are with updated parameters and with above two sub-plume scale corrections.

4. NCAR LES is shown just for reference. Flux profiles of ADC and Oceananigans correspond to grid scale only.