Possible addition of Brine plume parameterization

[vanroekel]

This potential feature is meant to mimic (in a simple way) the parameterization of Nguyen et al 2009 (https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2008JC005121) which is designed to limit the penetration of brine rejection plumes in models that often leads to overly deep haloclines in the arctic. It is my understanding of how Section 3.1 describes the parameterization. Translated to MPAS, I think this should be transporting the seaIceSalinityFlux which is the salinity rejected back to the ocean in ice formation. Here that flux is removed from the KPP non local flux and given a different shape function according to the paper. There is still work to be done, but was hoping for input on the implementation. In particular @maltrud, @milenaveneziani, @alicebarthel and @cbegeman I'd appreciate your thoughts.

[vanroekel]

My branch is farther along than ocean-discussions master so extra files are in the diffs. Here is a list of files I changed

mpas-ocean/src/Registry.xml
mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F
mpas-ocean/src/shared/mpas_ocn_tendency.F
mpas-ocean/src/shared/mpas_ocn_tracer_nonlocalflux.F

[vanroekel]

I've tested an initial implementation (currently in year 7), here is one plot of averaged MLD (density) in year 5 (JFM average)

There is a broad reduction in Arctic MLD.

[cbegeman]

@vanroekel I tested this in a column case and found that the salinity vertical mixing tendency is much more smooth through the boundary layer (in the presence of negative surface fluxes), which I think is what we want. (I think that there might be a difference in sign convention between the total tendency and the vertical mixing tendency in these figures.)

Brine parameterization off:

Brine parameterization on:

[cbegeman]

FYI, you may want to cherry-pick this commit to include the new config options in the build files: c83411f

[cbegeman]

Suggested change

	end if
	endif
	else
	tracersSurfaceFlux(index_salinity_flux, iCell) = tracersSurfaceFlux(index_salinity_flux, iCell) &
	+ seaIceSalinityFlux(iCell) * sflux_factor
	end if

I think we want to retain this (only relevant for diagnostic output?) because the surface flux is still seaIceSalinityFlux

[cbegeman]

Suggested change

	sis_flux(k) = min(1.0_RKIND,A*z**config_brine_param_n)*seaIceSalinityFlux(iCell)
	sis_flux(k) = max(0.0_RKIND,1.0_RKIND - A*z**config_brine_param_n)*-seaIceSalinityFlux(iCell)

When I scoped this out in python, a few other small changes were required to make your approach and mine equivalent (happy to share code), but the basic idea is to have sis_flux at the top of the model equal seaIceSalinityFlux and sis_flux below the BL equal 0. I think this is more intuitive because in its current form, it suggests that salt is fluxing from the bottom up rather than the top down, even though the divergence from both ways of writing it is identical.

[cbegeman]

Note: in Fortran 0**0 = 1 for most compilers so there would be an issue here if users set config_brine_param_n = 0. Probably the easiest thing to do is just set the value for z=0 outside this loop.

[cbegeman]

Just so I understand the relationship between our code and Nyugen, boundaryLayerDepth is determined based on the critical bulk Richardson number, right? Whereas Nyugen uses a threshold value of $d\rho/dz$.

[alicebarthel]

Is that the previous state boundaryLayerDepth? My main concern with this implementation is that if we are applying it after the OBL calculation in KKP, we are skipping one feedback Nyugen talked about: DeltaB increases --> OBL decreases.

[vanroekel]

my mistake I didn't realize the brine parameterization was detached from KPP. we definitely shouldn't use BLD here. In the paper is rho in situ density? Delta rho of 0.4 is quite large for potential density.

@alicebarthel can you explain how you would make the calculation on the current state? This call happens in the tendency part so it seems you'd want the previous (time level n) state for MLD. I feel i'm missing something.

[alicebarthel]

For readability, my preferred option would be to write the layer integral here. Any reason we prefer not to?

Suggested change

	A = (config_brine_param_n + 1) / D**(config_brine_param_n+1)
	k = minLevelCell(iCell)
	! need a salinity flux top and bottom and the flux into a cell is the divergence??
	do while (z < D)
	sis_flux(k) = min(1.0_RKIND,A*z**config_brine_param_n)*seaIceSalinityFlux(iCell)
	z = z + layerThickness(k,iCell)
	k = k + 1
	end do
	! need to pick up the last bit that goes over
	sis_flux(k) = 1.0_RKIND*seaIceSalinityFlux(iCell) ! this should finish
	do k=minLevelCell(iCell),maxLevelCell(iCell)
	! need to flip over so we take k+1 - k
	tend(index_salinity, k, iCell) = tend(index_salinity, k, iCell) + &
	sis_flux(k+1) - sis_flux(k)
	end do
	end if
	end do
	A = 1.0_RKIND / D**(config_brine_param_n+1)
	k = minLevelCell(iCell)
	! need a salinity flux top and bottom and the flux into a cell is the divergence??
	do while (z < D)
	sis_flux(k) = min(1.0_RKIND, A*((z + layerThickness(k,iCell))**(config_brine_param_n + 1) - z**(config_brine_param_n + 1)))*seaIceSalinityFlux(iCell)
	z = z + layerThickness(k,iCell)
	k = k + 1
	end do
	! need to pick up the last bit that goes over
	sis_flux(k) = min(1.0_RKIND, A*(D**(config_brine_param_n + 1) - z**(config_brine_param_n + 1)))*seaIceSalinityFlux(iCell)
	do k=minLevelCell(iCell),maxLevelCell(iCell)
	tend(index_salinity, k, iCell) = tend(index_salinity, k, iCell) + max(0.0_RKIND, sis_flux(k))
	end do
	end if
	end do

[alicebarthel]

Hi Luke, Carolyn,
A couple of thoughts on this.
i) I think writing the layer-integral form would be more readable than the flux divergence, unless I misunderstood something.
ii) is there a reason we prefer not to mimic the runoff spreading approach? It has the advantage of clearly being a surface flux spreading, rather than looking like double-counting (which I thought we were doing here at first glance).
iii) Does the use of boundaryLayerDepth means that we are skipping the KPP feedback that Nyugen was mentioning? It's still effective as is but I just want to clarify the choice that we make. Perhaps it's cleaner if the KPP code and plume parameterisation code remain independent?

[vanroekel]

Hi Luke, Carolyn, A couple of thoughts on this. i) I think writing the layer-integral form would be more readable than the flux divergence, unless I misunderstood something.

I don't have a strong opinion on this, either approach seems fine to me

ii) is there a reason we prefer not to mimic the runoff spreading approach? It has the advantage of clearly being a surface flux spreading, rather than looking like double-counting (which I thought we were doing here at first glance).

The advantage of this to me is it looks like the KPP non local flux which is a redistribution not double counting, when implemented correctly of course. But on the other hand it also seems pretty similar to spreading. Both are using some function to redistribute the surface flux. So I guess it doesn't matter to me which way to go and I think they are identical. If we move to spreading we need to correct the KPP surface buoyancy flux to account for the spreading. I don't think we want to use the same exponential spreading as is in runoff.

iii) Does the use of boundaryLayerDepth means that we are skipping the KPP feedback that Nyugen was mentioning? It's still effective as is but I just want to clarify the choice that we make. Perhaps it's cleaner if the KPP code and plume parameterisation code remain independent?

as I said I don't understand the concern. This (with fixes I realize now) would change the stratification even with the previous step BLD which would shallow the BLD in the next time KPP is called. I think it still hits that feedback. But I see use of BLD is not right anyway.

[vanroekel]

I think after quickly going back to Nguyen we should spread this over some MLD instead. Do you two think it's worth it to define a new MLD (threshold of 0.4) instead of just using the standard 0.03 de boyer montegut threshold?

[cbegeman]

B-case testing

in v3 SORRM configuration:
https://e3sm.atlassian.net/wiki/spaces/pd/pages/5655265281/20251107.v3.SORRME3r3.GMPAS-JRA1p5-DIB-PISMF.brineplume
run only 10 years because it enhanced MLD around with Weddell polynya

[cbegeman]

I think after quickly going back to Nguyen we should spread this over some MLD instead. Do you two think it's worth it to define a new MLD (threshold of 0.4) instead of just using the standard 0.03 de boyer montegut threshold?

Is the threshold of 0.4 the value of delta \rho from the surface density? Going back to the Duffy et al 1999 approach cited in the paper?

[vanroekel]

I think after quickly going back to Nguyen we should spread this over some MLD instead. Do you two think it's worth it to define a new MLD (threshold of 0.4) instead of just using the standard 0.03 de boyer montegut threshold?

Is the threshold of 0.4 the value of delta \rho from the surface density? Going back to the Duffy et al 1999 approach cited in the paper?

Yes sorry I didn't realize initially I cited the Duffy 1999 (and is relative to the surface density) approach, I see the d rho/dz criterion in the Nguyen paper as you mentioned. They use a threshold of 0.005 (a couple of paragraphs into section 3.1). We do have this calculation in the mixed layer AM so could use that.

[cbegeman]

I think after quickly going back to Nguyen we should spread this over some MLD instead. Do you two think it's worth it to define a new MLD (threshold of 0.4) instead of just using the standard 0.03 de boyer montegut threshold?

Is the threshold of 0.4 the value of delta \rho from the surface density? Going back to the Duffy et al 1999 approach cited in the paper?

Yes sorry I didn't realize initially I cited the Duffy 1999 (and is relative to the surface density) approach, I see the d rho/dz criterion in the Nguyen paper as you mentioned. They use a threshold of 0.005 (a couple of paragraphs into section 3.1). We do have this calculation in the mixed layer AM so could use that.

I don't have any insight into which approach is better off hand. If it's not too noisy, the drho/dz threshold seems more physical.

[vanroekel]

I don't have any insight into which approach is better off hand. If it's not too noisy, the drho/dz threshold seems more physical.

I think I'd lean to gradient too if we are looking for a pycnocline as the paper suggests. Maybe I can turn that gradient MLD calculation on in my current spin up run to see how it compares with the threshold in particular if it is noisy.